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SUMMARY 

A thin film of oil on a surface responds primarily to the wall shear stress generated 
on that surface by a three-dimensional flow. The oil film is also subject to wall pressure 
gradients, surface tension effects and gravity. The partial differential equation governing 
the oil film flow is shown to be related to Burgers 1 equation. Analytical and numerical 
methods for solving the thin oil film equation are presented. A direct numerical solver is 
developed where the wall shear stress variation on the surface is known and which solves 
for the oil film thickness spatial and time variation on the surface. An inverse numerical 
solver is also developed where the oil film thickness spatial variation over the surface at 
two discrete times is known and which solves for the wall shear stress variation over the 
test surface. A One-Time-Level inverse solver is also demonstrated. The inverse numerical 
solver provides a mathematically rigorous basis for an improved form of a wall shear stress 
instrument suitable for application to complex three-dimensional flows. To demonstrate the 
complexity of flows for which these oil film methods are now suitable, extensive examination 
is accomplished for these analytical and numerical methods as applied to a thin oil film in 
the vicinity of a three-dimensional saddle of separation. 
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1. INTRODUCTION 


The response of a thin oil film on a surface subjected to a three-dimensional (3D) 
aerodynamic flow proves to be of interest to the fluid mechanics community for several 
reasons- First, the study of thin oil films has led to new instrumentation which accurately 
measures the wall shear stresses generated by 3D flows. Second, extensive use of oil flow 
visualization on surfaces bounding 3D flows requires that we understand the limitations 
of these oil flow techniques in the neighborhood of the complex topological features, in- 
cluding singularities, of the limiting surface streamlines. Additionally, the study of thin 
oil films provides an opportunity to consider the solution of nonlinear partial differential 
equations with clear practical significance and which require minimal computational and 
programming resources. 

Instrumentation based on thin oil films has undergone sustained development (Tanner, 
et al. 1-4 ; Monson, et al. 5-7 ; Mateer, Monson, and Menter 8 ) to provide accurate measure- 
ment of the wall shear stresses generated on a test surface by aerodynamic flows. Recent 
improvements (Naughton and Brown 9-11 ), particularly for 3D flow applications, in the 
form of the oil film wall shear stress instrument have been made possible by the appli- 
cation of Computational Fluid Dynamics (CFD) solution techniques to the thin oil film 
equation. Experimental wall shear stress measurements guide development of turbulence 
models to further improve CFD solvers. An important barrier to the accurate prediction 
by Navier-Stokes CFD solvers of complex 3D flows is the accurate modeling of the flow- 
field turbulent Reynolds stresses. Considerable effort has been expended over the past six 
decades in turbulence modeling with advances in this difficult area being frustratingly slow. 
The pace of improvement in turbulence modeling has improved recently, however, partly 
due to a maturation in the numerical methods of Navier-Stokes solvers, and partly due to 
enhanced instrumentation. Two instruments most relevant to improvements in turbulence 
modeling are the laser Doppler velocimeter (LDV) and the laser interferometric skin fric- 
tion (LISF) instruments. The LISF and related successor skin friction instruments deduce 
the wall shear stress by analysis of either the thinning rate or the detailed shape of an oil 
film which is initially spread on a test surface and which then responds to the aerodynamic 
flow over that test surface. Prior to this study, self-similar ID solutions for the oil film 
response were used 1-8 in the analysis for the oil- film based instruments. However, as these 
oil-film based instruments are now being applied to complex 3D flow situations, a more 
rigorous treatment of the oil film response is now required and is addressed by this present 
work. 

A further reason to study the thin oil flow equation is the utility of surface oil flow vi- 
sualization technique (Maltby 12 ). Surface oil flow visualization is one of several techniques 
to discern the limiting streamline flow patterns on the surface. An aid to understand- 
ing complex 3D flows is the topological analysis of these limiting streamlines (Tobak and 
Peake 13 ; Chapman and Yates 14 ). Visualization, both of numerical solutions and experi- 
mental flows, of these surface streamlines and associated singularities provides an overall 



topological framework for categorizing and comprehending the complex flow patterns which 
may arise. For 3D flows, the types of singularity points which may occur on surfaces are 
the saddle of separation, the saddle of attachment, the node of separation, the node of 
attachment, the focus of separation, and the focus of attachment. Limiting separation or 
attachment lines connect these point singularities. A further surface streamline topolog- 
ical rule is that these singularity points appear in combinations such that for a simply 
connected closed body: 

N-S = 2 

where N is the total number of nodal and focal points appearing in the flow and S is 
the total number of saddle points appearing in the flow. Easily overlooked are the node 
of attachment at the nose of an aerodynamic body and the node of separation at the 
tail. Symmetry of the flow can also lead to an undercount since a singular point may then 
actually appear twice. The flowfield need not be uniquely defined by the surface streamline 
patterns observed. 

The importance of the ability to correctly identify these singularity points and further 
to predict with a Navier-Stokes solver the location of these singularity points is generally 
underestimated. In particular, the saddle of separation and the saddle of attachment can 
be difficult to identify and yet either can appear in some flows at a particular location with 
significant impact on the flowfield pattern above the surface. To miscalculate a singularity 
point for a turbulent flow may be the consequence of an improper grid or more importantly 
an improperly constructed turbulence model. 

Squire 15 provided the first theoretical study of the thin oil film equation. Numerical 
solutions for the ID thin oil film under a 2D aerodynamic flow were presented. In partic- 
ular, the behavior of oil flow visualization in the vicinity of 2D separation was addressed. 
Squire concluded that the thin oil film did not significantly alter the boundary layer and 
that oil flow visualization would tend to indicate 2D separation slightly upstream of the 
actual boundary layer separation due to pressure gradient effects. Part of the purpose of 
this present work is to extend the earlier work of Squire to consider 2D thin oil films under 
3D aerodynamic flows. 

A third reason for the study of the thin oil film equation lies in its utility in the study 
of numerical methods. The thin oil film equation is an extension of the familiar Burgers 
equation often used to test CFD numerical methods. As with Burgers’ equation 16 , the thin 
oil film equation is a scalar hyperbolic wave equation which may be solved by numerous 
solution methods, including finite-difference, finite- volume, characteristic and Lagrangian 
techniques. These thin oil film solutions can be accomplished in ID or 2D on a workstation. 
Additional advantages accrue to the study of the thin oil film equation, how'ever, in that 
for both the ID and 2D forms the eigenvalue(s) and hence characteristic direction can be 
forced to change sign at a particular location in a model problem. Further, in the present 
article, we will demonstrate an inverse solution numerical method useful for the skin friction 
instrument. Additionally, the thin oil film equation offers a context to study ID and 2D 
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model problems both numerically and experimentally which can provide particular insight 
when fluid mechanics students are introduced to numerical methods. 

Note that the ID thin oil film problems are associated with 2D aerodynamic flows, 
while 2D thin oil film problems are associated with 3D aerodynamic flows. Where clarity 
in this paper requires, we will specifically state, for example, 2D thin oil problem or 
2D aerodynamic flow. Also, the two fluids need not be restricted to oil and air, but to 
clearly distinguish between the flow of the two fluids, we shall use the terms “oil flow” and 
“aerodynamic flow” throughout the remainder of the paper. 

In the present study, both direct and inverse numerical solution techniques for the 
thin oil film equation are developed. The direct numerical solver considers the case where 
the wall shear stress field on the surface, J^ w (x, 2 ), is known and the direct solver provides 
the oil film thickness variation with time over the surface, h(x,z,t). The inverse numer- 
ical solver considers the case where the oil film thickness variation at two discrete times, 
h{x , 2 , ti) and h(x, z , £ 2 ), is known and the inverse solver then provides the wall shear stress 
variation over the surface, ~r* w (x. z). 

In the sections to follow, the thin oil film equation is first derived. Next, various 
solutions are demonstrated including exact ID and 2D self-similar cases. Then numerical 
procedures for both the direct and inverse solutions are presented. These solvers are then 
applied to several example practical problems. Programs to solve these example problems 
are written in the c programming language for use on a Unix workstation and are available 
from the first author. 
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2. THE GOVERNING PARTIAL DIFFERENTIAL EQUATION 


The Navier-Stokes equations describe the response of a thin film of viscous liquid, 
typically oil, which is initially spread on a surface and which then experiences a 3D flow 
of a second fluid, typically air, over that surface. Simplifications to the Navier-Stokes 
equations for such an oil film flow are possible which still result in an accurate description 
of the oil flow while considerably reducing the difficulty of the solution method. The thin 
oil film equation derived below is essentially the continuity equation integrated across the 
thickness of the thin oil film, with additional information incorporated from simplifications 
of the x- and z-momentum equations. The derivation of the thin oil film equation is quite 
straightforward but is described below to establish the restrictions on the equation and to 
clarify the numerical procedures used to solve the equation. 

Consider a thin film of viscous liquid, such as a silicone oil, initially placed on a test 
surface as shown in figure 1. Typically, the thickness of oil is a few microns in thickness 
which varies with location and time. This test surface, and thus the oil film, is then 
subjected to a 3D aerodynamic flow over the test surface. The 3D aerodynamic flow 
generates on the test surface a wall shear stress vector, 7 * = (t x (x, z), t z (x, z)), acting 
tangential to the surface, and a wall normal stress or pressure, P(x , z), acting normal to 
the surface. The oil film will flow in response to these wall stresses and to the gravitational 
body force acting on the oil. Additionally, the oil film will experience surface tension 
effects related to the curvature of the oil film surface. For the purposes of this derivation 
we assume the 3D aerodynamic flow is steady with time. 

The thickness of the film, h(x,z,t), will vary with position on the surface and with 
time. To derive the differential equation governing the oil film behavior, consider the 
control volume of figure 1. The control volume encloses the full height of the oil film, h, 
and is of finite length, Ax, and width, Az, aligned with the x and z axes, respectively. 
Thus, the oil mass in the control volume at any time is given by m cv = p a hAxAz. A 
change in film thickness, h , and, thus, mass in the control volume occurs during a time 
interval, A due to the differences in mass flux normal to the four sides of the control 
volume through which oil may flow: 

Am cv / At -I - A X F T A Z G — 0 


where 


m cv — pohAxAz 

-h 


= J j p 0 udA = J p 0 udyAz = p 0 U c h,Az 
— J j p 0 wdA = J p 0 wdyAx = p 0 W c hAx 
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Taking limits of Ax, A z and At => 0, we obtain: 


dh/dt + dU c hjdx + dW c h/dz = 0 


where we have defined: 


U c h = I udy , and W c h = I wdy 

Jo Jo 

The fluxes, F and G, can be evaluated by means of a low-Reynolds number simplifi- 
cation of the x- and z-momentum equations. The film Reynolds number may be evaluated 
as a ratio of inertial effects to viscous effects: 

Re f = (p 0 U 2 /L)/(p 0 U c /h 2 ) = (p 0 U c h/p 0 )(h/L) 

To estimate the film Reynolds number we make use of Tanner’s ID self-similar solution 
result, U c = rhf2p 0 , along with estimates for h = 10 _6 ra, L = 10 ~ 2 m, r = 20 N/m 2 , 
p 0 = 1000 Kg/m 3 , and u Q = 100 centiStokes to obtain: 

Re/ ~ rh z j2p 0 u 2 L = 10 -8 

As a consequence of the low Reynolds number, we ignore the inertial terms in the x- 
momentum equations, giving within the oil film: 


0 = dr Xf0 /dy - dP Q /dx + p Q g x 

For the purposes of clarity in this derivation, we introduce the subscript, o, to signify the 
shear stress and pressure within the oil film, and the subscript, a, to signify the aerodynamic 
wall shear stress and wall pressure which are applied as boundary conditions to the oil film 
at the air/film interface located at y = h. 

Integrating, from the air/film interface inward, for the shear stress variation through 
the film layer: 



~< dP o 

Tx.o T x,a — v q 


Po9x){y h) 


Note, the y-momentum equation implies that the pressure. P 0 (x, z), may be assumed 
constant across the oil film thickness. 

Integrating again, but from the wall out into the film, gives: 


/ y rv du dP a 

Tx.ody = y o = VoU = T x ,aV + - Po9x)(y 2 /2 - hy) 

Or, considering both the x and z components of velocity within the oil film: 
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(2.1a) 


u = (T x , a y + - Po 9 x){y 2 /2 ~ hy))/p 0 

ftp 

W = (r z , a y + - p 0 gz)(y 2 / 2 - hy))/p Q (2.1fc) 

Integrating yet again finally gives us an expression for the mean convective velocity: 

rh Qp 

U c h = J udy = T x , a h 2 / 2 p 0 - (-^7 - Po 9 x){h z / 2 >p 0 ) 

Similarly, from the 2-momentum equation: 

rh Qp 

W c h = J wdy = T z , a h 2 / 2 p 0 -{-^-po 9 z)(h 3 / 3 p, 0 ) 

To account for surface tension effects, note that the pressure, P 0 , within the film will be 
altered from the aerodynamic wall pressure, P a , due to the curvature (1 /R x and 1/P Z ) of 
the oil film surface, hence: 


P 0 — P a + <7(1/ R x + 1/ R z ) ~ P a &{h X x d" ^zz) 

Summarizing, the differential equation governing the response of a thin film of oil to 
a 3D aerodynamic flow is: 



W c = 


r z h 

2 p 0 


dh dUJi dWgh 
dt dx + dz 

, dP _ da{hxx + hzz ) 
dx dx 


, dP _ dajhxx + h Z2 ) 
dz dz 


= 0 

Po9x){h /3/i 0 ) 
Po 9 z)(h /3/x 0 ) 


(2.2a) 
(2.2 b) 


(2.2 c) 


Equation 2.2, with slight rearrangement, was first given by Squire 15 , and we refer to this 
equation as “Squire’s Form” of the thin oil film equation. Tanner 1 gave a different form 
which we refer to as “Tanner’s Form” of the thin oil film equation. In Appendix A, the 
equivalence of the two forms of the thin oil film equation is demonstrated through a metric 
transformation. 

The issue of boundary conditions will be treated in the solution subsections below. 
Note, we have dropped the subscript, a. on r x , r z and P in the equation above and for the 
rest of the paper since this subscript was introduced for clarity to distinguish between the 
values within the oil film and the values imposed from the aerodynamic flow. Henceforth, 
the aerodynamic wall shear stress and wall pressure meanings for these terms are assumed. 
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In the absence of the surface tension term, which for many problems is negligible, the 
thin oil film equation is first-order hyperbolic or wave-like with the characteristic direc- 
tion of information propagation being indicated by the mean convective velocity vector, 
( U C1 W C ). With inclusion of the surface tension term, the equation becomes elliptic. 

Further observe that the r terms are multiplied by h 2 whereas the remaining terms 
are multiplied by h 3 . Thus, for h very small, the r terms typically dominate except near 
singularities, where the shear stress approaches zero. 

For the case where the pressure, gravity and surface tension terms are negligible 
compared to the shear stress term, the thin oil film equation becomes: 


dh drJfJ2±o , chj}fJ2±o _ n 

dt dx dz 

In a coordinate independent form, the above equation becomes: 


( 2 . 3 ) 


dh . h 2 

dt 2 Ho 


7 *) = 0 


( 2 . 4 ) 
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3. ANALYTICAL SOLUTIONS 


Analytical solutions for the thin oil film equation are considered in this section. These 
analytical solutions are for relatively simple test case conditions, but do include both ID 
and 2D thin oil film problems for 2D and 3D aerodynamic flows, respectively. These 
analytical solutions can be instructive as to general thin oil film behavior as well as being 
useful in their own right. More general cases, whether for ID or 2D oil films, require 
numerical solution procedures. Even for numerical procedures, the analytical solutions aid 
in formulating boundary conditions. Furthermore, the analytical solutions considered in 
this section then provide known test cases to assess the validity and accuracy of the more 
general numerical solution procedures. 

We start by establishing a self-similar form of the 2D oil film equation. Then, we 
further reduce this equation to an ordinary differential equation for the ID self-similar thin 
oil film problem, and establish several self-similar relationships. Analytical ID solutions 
are then established. We then return to the 2D self-similar form, and, for special forms of 
the wall shear stress under a 3D aerodynamic flow, establish 2D analytical solutions. 

Consider the thin oil film equation where the pressure gradient, gravity and surface 
tension terms are negligible: 


dh d ,T x h 2 d r z h 2 

dt + dx 2fi 0 dz 2 Ho 


) = 0 


We consider solutions which are self-similar in time of the form: 


(3.1) 


h(x,z,t) = H(x,z)/t (3.2) 

Thus, dh/dt = - H/t 2 , giving the 2D self-similar form of the thin oil film equation: 


_ s +°&!L )+ °e£. ) = o 

ox 2 fi 0 oz 2 Ho 

In coordinate independent form, the above equation becomes: 


(3.3) 


-ff+V.(^7») = 0 (3.4) 

The self-similar solution, H{x,z), describes the asymptotic shape of the thin oil film at 
large time. 


3.1 ID ANALYTICAL SOLUTIONS 

For a ID thin oil film, the self-similar partial differential equation reduces to an ordi- 
nary differential equation: 
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One ID self-similar relation (Tanner and Blows 1 ) can be established by integrating 
the above equation from the leading edge where H = 0 at x = xo: 

i£ Bdx -m£ h *‘ < 3 ' 6) 

Another self-similar relation is found by defining ( = ( r/fi 0 ) l / 2 H and substituting in 
equation 3.5 giving: 

C Wr) ,/! -^K !/ 2) = ° 

or, rearranging 


0u o /r) 1/2 = dC,/dx 

Integrating, with Q = 0 at x = xo, and rearranging gives the ID self-similar relation (also, 
Tanner and Blows 1 ): 


H = ht = (h 0 /t) 1/2 f (fi 0 /T) l/2 dx (3.7) 

Jx 0 

The relation given by equation 3.7 may be solved for H(x) given known t{x) by 
numerical methods or, where suitable, in closed analytical form. Likewise, the relation 
given by equation 3.6 may be solved for r(x) given known H(x). The studies of Tanner 
and of Squire provide several such ID solutions. 

The simplest ID self-similar thin oil film solution is for the case of constant wall shear 
stress, t(x), where: 


Hr = htr = fj, 0 (x — xq) 


(3.8) 


Axisymmetric Analytical Solutions 

For flows over axisymmetric bodies, the governing equation for a thin oil film becomes: 

dh l c? rrh 2 

m +( -r ) rs^ ) = 0 (3 - 9) 


and the time self-similar form becomes: 



(3.10) 


H-{hU ^)-0 

K r ds K 2 fi 0 

Rearranging and integrating from sq, where h = 0, to s: 


T =^f K rHds= ^[j his (3I1) 

Alternatively, substituting £ 2 = ttH 2 j \i Q into equation 3.10: 

H = ht = ds (3.12) 

rT Jso T 

As one example of a closed form axisymmetric thin oil film solution, consider the case 
where the aerodynamic flow consists of the region about an axisymmetric stagnation point 
or node of attachment formed by placing a circular plate normal to the flow. For this 
case, r = s. A laminar solution for this axisymmetric stagnation point has been given by 
Homann 17 , and is reported in both White 18 and Churchill 19 . The wall shear stress on the 
plate varies linearly as: 


r = /3r, (3>0 (3.13) 

Assume the “leading edge" of the axisymmetric oil film is located some small distance, 
r 0 , from the stagnation point. Integrating equation 3.12 with the known wall shear stress, 
equation 3.13, gives the oil film shape: 

H — ht — (/z o /0)( 1 - r 0 /r), r > r 0 (3.14) 

For distances far from the stagnation point the oil will tend toward a uniform thickness, 
which varies inversely with time. For locations close to the stagnation point, pressure 
gradient effects become important which then require a non-similar numerical solution. 
For an axisymmetric node of separation, we may also assume that the local shear stress 
varies according to r = /3r, except (3 < 0. The oil film leading edge is applied at r — 7~o, 
and the flow of oil is inward toward the node of separation. The solution then is: 

H = ht = (vo/miro/r - 1), r < r 0 (3.15) 

For the axisymmetric node of separation, the singularity at r = 0 is avoided due to 
both surface tension and pressure gradient effects not included in this time self-similar 
solution. 
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3.2 2D ANALYTICAL SOLUTIONS 


Closed form analytical solutions may also be found for a number of interesting cases 
of a 2D thin oil film responding to a 3D aerodynamic flow. The flow of a 2D thin oil film 
in the immediate vicinity of 3D surface streamline topological singularities, such as nodes 
and saddles, leads to closed form solutions. Additional 2D thin oil film test cases suitable 
for testing numerical 2D thin oil film solvers involve axisymmetry and are also treated in 
this section. These closed form solutions serve as suitable test cases for the more general 
numerical 2D thin oil solvers, both direct and inverse, discussed in following sections. In 
this subsection, we develop several of these closed form analytical solutions. 

Analytical Solutions for Saddles and Nodes 

A linearized form of the flow field about surface streamline singularities in 3D aerody- 
namic flows was developed by Perry and Fairlie 20 . In this section, we make use of relations 
based on their work to provide appropriate surface shear and wall pressure fields which 
may be generated by the fiowfield so as to study the response of a thin oil film about 
such surface streamline singularity points. Through a suitable coordinate rotation and 
stretching 20 , a “canonical” form of the surface singularities may be obtained. Thus, in 
this section, although we consider saddles located on a plane of symmetry, the results are 
more general. The papers of Perry and Fairlie 20 and of Hung, Sung, and Chen 21 should 
be referred to for a more extensive development of the shear field and pressure field in the 
vicinity about a singularity or critical point in 3D aerodynamic flows. 

Consider the case of a surface streamline singularity, either a saddle or a node, located 
on a surface along a plane of symmetry. Figures 2a and 2b show the defining streamlines 
about a saddle of attachment and a saddle of separation, respectively. The streamlines 
depicted are either constrained to the surface or form in a symmetry plane for the aero- 
dynamic flow above the saddle. Note for the saddle of separation, a streamline originates 
from the singularity point and departs into the flow upward from the surface. For a saddle 
of attachment this streamline will have flow toward the surface. Figures 3a and 3b further 
depict the surface streamlines about a saddle of attachment and a saddle of separation, 
respectively. Figures 3c and 3d depict the surface streamlines about a node of attachment 
and a node of separation, respectively. For convenience, we place the origin of the coordi- 
nate system at the location of the surface streamline singularity, ( x,z ) = (0,0), with the 
plane of symmetry occurring along 2 = 0. At the surface singularity point, both the x- 
and 2 -component of the wall shear stress go to zero, ((r x ,r z ) = (0,0) at (x.z) = (0,0)). 

For some region locally about this singularity point, the wall shear stress varies ac- 
cording to: 


'T x — (%3C 
t z = bz 


(3.16) 


11 



The relative values of a and b determine the characteristics of the singularity. The 
saddles and nodes of separation and attachment may thus be distinguished by the values 
of a and b : 


s <C 0, 6 < 0. a + b <C 
a < 0, 6>0, a T b < 
a<0, 6 > 0, a + b > 
a>0, 6 > 0, a + b > 


0 : node of separation 
0 : saddle of separation 
0 : saddle of attachment 
0 : node of attachment 


(3.17) 


The pressure field on the surface may also be found, following the analysis of Perry 
and Fairlie 20 : 


dP 

— — = —(3a + 6)/ tan# 
ox 

f- 


(3.18) 


where 9 is the departure angle of that limiting streamline which emanates from the singu- 
larity point and which departs the surface into the flow above as depicted in figure 2b. As 
an aside, equation 3.18 implies that for a given value of dr x /dx = a < 0, a lower pressure 
gradient will be associated with a saddle of attachment than the pressure gradient associ- 
ated with a saddle of separation. This may be of practical importance where a laminar flow 
generates a saddle of attachment for a given geometry, while a turbulent flow generates a 
saddle of separation for the same flow geometry. 

For the two types of saddles and for the node of separation, we have a < 0. The oil 
leading edge is applied as a straight edge at, for example, x$ < 0, and the oil flows toward 
the saddle point or node. The solution domain for the thin oil film is then xq < x < 0. 
For the node of attachment, flow is away from the node, and we may chose the solution 
domain as x < x 0 < 0. For the node of attachment, the leading edge of the oil film can 
also be chosen to be at the node point (xo — 0)- Note that the special case of either the 
axisymmetric node of separation or of attachment as treated in the previous section differs 
in that a = b and the oil leading edge is applied not as a straight line but at the circle 
defined by r = ro- 

We first derive the shape of the surface streamlines, since along any given surface 
streamline: 


dz t z _ bz 

dx t x ax 

If we assume the streamline passes through the arbitrary point, (x s , 2 S ), we may integrate 
the above streamline relation to obtain: 


12 



z 

Zs 


(_£.)&/<* _ y>/-a 

X 


(3.19) 


Curves demonstrating the above streamline equation for several values of a/b are shown 
in figures 3a-3d for the saddles and nodes of attachment and separation. 

To proceed toward a closed form solution for the thickness distribution for an oil film 
near a surface streamline singularity point, first consider the saddle of separation. The 
streamline flow along the surface in the plane of symmetry is toward the saddle, and thus: 


a < 0, and b > 0 

Consider that the wall shear stress effects predominate over the pressure gradient and 
surface tension effects, except very close to x = 0. To more properly examine the inclusion 
of these effects, we will later resort to numerical solvers. With the assumption of negligible 
pressure and surface tension, the 2D thin oil film equation takes the form of equation 3.1, 
which we repeat here: 


dh d ,r x h 2 d ,r z h 2 
Tt + Tx ( ^: ) + Tz ( = 0 

First, apply the known shear field of r x = ax, and r z = bz , and expand: 


dh 

dt 


d xh 2 

dx 2/i 0 * dz x 2fx 0 


d .h 2 , L/ 
z~(7r~) + K^~) = o 


2 Ho' 


Note that if the oil thickness, h, initially varies at most in x then there is no driving 
force for h to subsequently acquire a variation in the z-direction. Thus, dhjdz = 0 not 
only initially, but for all time. A variation in x must occur due to the leading edge. Thus, 
h = h(x,t), which leads to: 


dh d xh 2 h 2 

dt + a dx 2 Ho + 2 Ho 


) = 0 


Now, applying the self-similarity relation (h(x,t) = H(x)/t) and rearranging (b ^ 0 
avoids the trivial ID case): 


dH/(ft 0 — (a + b)H/2) = dx/(ax) (3.20) 

Excluding the special case where a + b = 0, we integrate from xo, where h = 0. to x: 

H = (— ^r)[l - ( — )'- 2q) ], for 0 < x < x 0 (3.21) 

G i 0 X q 

This solution, also see Tanner 4 , covers the interesting case of thin oil film response 
in the vicinity of the surface singularities known as saddles of separation and saddles of 
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attachment and is also valid for the node of separation, when such singularities are located 
on a plane of symmetry. We will in later sections of this paper also obtain numerical 
solutions and discuss in greater depth. 

A node of attachment (a + b > 0, a > 0, b > 0) has a different domain of solution and 
thus different limits of integration for equation 3.20, resulting in: 

H = — t)[ 1 - for 0 < x 0 < x (3.22) 

cl H- b x 

Oddly, for a node of attachment with the “leading edge” located at xo = 0, the only 
allowable self-similar solution is: 

H = 2 Ho/ (a + 6), for all x 

The behavior of the oil film at the location of the saddle differs considerably for the 
two types of saddles: 


saddle of separation(a + f> < 0, a<0, 6 > 0) : H — > oo, as x — » 0 
saddle of attachment^ + b > 0, a<0, 6 > 0) : H = 2n 0 /(a + 6), for x = 0 


Al$o, 

node of separation(a -I- 6 < 0, a < 0, 6 < 0) : H — » oo, as x — * 0 

For the special case where a + b = 0, note that for oil to flow from the oil film leading 
edge at x = xo to the singularity at x = 0, we must have a < 0. For this special case, 
integration of the ODE above from xo toward x = 0 yields: 

ht = H = (-Ho/ a) lnxo/x, for x < xo, a < 0, and a + b — 0 (3.23) 

Examination of the form of these solutions, equations 3.21 and 3.23, suggests that a 
suitable plot of H vs logio x should prove useful in determining the ratio of the wall shear 
stress slopes b/a and thereby, for example, allow determining whether a saddle is a saddle 
of attachment or saddle of separation. 


Analytical 2D Thin Oil Film Relations 


For somewhat more general 2D thin oil film cases we can establish useful relations. 
Consider equation 3.1, where we further assume that t x is a function of x only, and 
t z = bz + /(x): 


dh 

dt 


d ,r x h 2 d h 2 h 2 

+ 7T(t— + r *7T T - T“) 

dx 2 Ho oz 2/i 0 2 ii Q 


(3.24) 


14 



Note that h = h(x,t) is a possible solution. If the initial h is only a function of x , 
(h 0 = h 0 (a:)), then initially d(h 2 )/dz = 0 and, consequentially, there is no driving force 
for h to subsequently acquire a variation in the 2 -direction. Thus, we may seek solutions 
for h where h is only a function of time and x. 

Applying the self-similarity relation where h = H/t, noting d/dz = 0. and rearranging: 


d ,t x H 2 

— (— — 
dx 2 Ho 


) = H 


Integrating from xq, where h = 0, to x: 


bJP 

2 ^ 


(3.25) 


T x = 



< hH ^ 

( Wo )]dX 


(3.26) 
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4. NUMERICAL SOLUTIONS 


Numerical procedures are required to solve for the general time-dependent response 
of a thin oil film on a surface subjected to an aerodynamic flow. In the previous section, 
we discussed special cases where time self-similar solutions of analytical form were pos- 
sible. However, the inclusion of initial conditions, pressure gradient, gravity, centrifugal, 
or surface tension effects, lead to a time-dependent oil film response which needs not be 
self-similar in time. Also, the geometry of the oil film or the aerodynamic flow will not 
typically lend itself to a known analytical solution. Fortunately, CFD methods are readily 
available for application to the numerical solution of the general thin oil film problem. In 
this section, we describe those methods which we have found most suitable. 

Two primary numerical approaches were developed in the present studies. The first 
approach we refer to as a direct solver, where the oil film thickness, h(x, z , t), is solved for 
knowing the aerodynamic wall shear stress and wall pressure, and where gravity, centrifugal 
and/or surface tension effects are included. The second approach we refer to as an inverse 
solver, where the wall shear stress, (x, z ), is deduced knowing the response of the oil film 

thickness at several times, h(x, z , fi) and h(x, 2, £2)1 and the surface flow direction, 7(2:, 2), 
as well as known wall pressure, and where the gravity, centrifugal and/or surface tension 
effects are included. The inverse solver, in particular, provides a rigorous foundation for 
the oil film method of experimental measurement of wall shear stress. 

The numerical method described is chosen for its simplicity for use with both the 
direct and inverse solvers. For the direct solver, the h solution is advanced in time, with 
repeated sweeps through the grid. For the inverse solver, the r solution is simply marched 
once in space through the grid. The dominant physical phenomenon described by the thin 
oil film equation is hyperbolic, having a known characteristic direction. The characteristic 
direction may be thought of as the direction in which information propagates. Each point 
in the thin oil film is only under the immediate influence of its upwind neighboring points. 
Thus, the dominant hyperbolic nature of the thin oil film equation allows the direct solver 
to advance the solution in time by means of a point-by-point implicit solution. Each point 
is advanced in time once those neighbor points which are upwind are advanced in time. 
Each time-advance sweep through the grid proceeds from the boundary points inward to 
the interior points. 

For the direct solver, addition of the surface tension term introduces an elliptic feature, 
which allows information to propagate in all directions. However, the influence of the 
surface tension effects on a thin oil film will occur over only a quite limited region. The 
oil film is thin (typically a few microns) and is assumed to cover an extended region of the 
test surface. For surface tension to be significant, the curvature of the film surface must 
be significant. For the film to remain thin, this curvature cannot cover a large region. A 
variety of elliptic solution methods may be applied to implement the surface tension terms. 
Particularly in two directions, implicit methods unnecessarily complicate the solver. Due 
to the large time scale associated with the surface tension terms, however, the simplest 
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approach is to incorporate the surface tension explicitly and limit the time step. The time 
scale for surface tension, (A t a % fji 0 L*/ah 3 ), is typically quite large (days rather than 
seconds), particularly when compared to the time scale required for accurate solution of 
the shear stress terms. Thus, explicit treatment of the complete surface tension terms in 
an otherwise implicit point-wise direct solver should prove sufficient and quite practical 
for most shear-driven applications. An example of a line-implicit treatment of the surface 
tension term (for one predominant direction only) is included in the application section. 

For the inverse solver, in contrast, the thin oil film equation remains hyperbolic, even 
with the inclusion of the surface tension terms, which then are known source terms. 

In the following subsections, we first describe the numerical procedures we use for an 
interior node of the direct solver. A Box-Implicit numerical method is described. The as- 
sociated numerical boundary conditions axe then described. An alternative Finite- Volume 
Upwind-Implicit numerical method is also presented. In the final subsection, the inverse 
solver numerical methods axe described. 

4.1 DIRECT NUMERICAL SOLUTIONS 

To simplify the derivation of the numerical method, we first neglect the pressure 
gradient, gravity, centrifugal and surface tension effects and then add these effects later 
in the section. Considering only the viscous terms in the oil flow, the partial differential 
equation for an interior point of a thin oil film is given by a simplified version of equation 
2 . 2 : 


dh 

dt 


+ l (CW + l (wy, >=° 


(4.1a) 


U c = 


Tx h 

2 Vo 


(4.16) 


W c = 


Tjh 
2 fi 0 


(4.1c) 


The above partial differential equation is hyperbolic and at any point considered 
the equation represents a specialized form of the continuity equation which depicts the 
convection of the conserved variable h (related to the oil mass at a point) at a velocity 
U c in the xr-direction and at a velocity W c in the z-direction. Any numerical method to 
be successful for application with this equation should be conservative and time-accurate, 
with consideration given to the direction of the convective velocity. 

First integrate the partial differential equation at an arbitrary point with respect to 
time, between the two times, t\ and 1 2 , giving: 




(4.2) 
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Note, the wall shear stress components, r x (x, z) and r z (x, z), are assumed to be known 
from the aerodynamic flow and are steady in time. We also assume the oil viscosity is 
steady in time. We have seen in the prior section on analytical solutions that the oil film 
thickness tends to vary inversely with time according to h ~ t~ l . Because of this inverse 
time relation over much of the oil film, greater practical accuracy can be obtained by 
evaluating the integral with the numerical quadrature given by: 

f h? dt — hih,2(t2 — £ 1 ) (4-3) 

Jt i 

The above quadrature is exact for h « £ _1 while only formally first-order accurate. 
Use of a formally second-order accurate quadrature, such as (hi + fi 2) 2 (£2 — U)/4, actually 
results in a reduction in time accuracy for many practical cases where the h « £ -1 time- 
similarity is approached. As a consequence of the time-integration, equation 4.2 becomes: 

(/*2 — h\) + —(—2-h\h2&t) 4- — (— hih 2 M) (4-4) 

dx 2 Ho dz 2/z 0 


Box-Implicit Direct Solver: Interior Node 

We now require spatial-discretization of the oil film into a 2D array of nodes, (i = 
l,zmax, and j = 1 , j max) , with even spacing, Ax and A z in x and z, respectively. Thus, 
Xij = iAx and Zij = jAz. More complex gridding treatments can be accommodated but 
are not treated here. 

To develop the Box-Implicit numerical procedure, consider the solution molecule 
shown in figure 4. It is assumed that we know the solution at the time level <1 and 
that we desire to advance the solution at node (z, j) to the new time level, £ 2 - We assume 
the convective velocities are U c > 0 and W c > 0. Thus, to solve for /&2,i,j = h(xi^, Zij.t 2 ) 
we require the solution at nodes (i — l,y), (t, j — 1), and (i — 1, j — 1) to have already been 
advanced to the new time level, £ 2 - To achieve second-order accurate approximations to 
the spatial derivatives, d/dx and d /dz, we evaluate at the midpoint of the control volume, 
i — 1/2, j — 1/2: 


dx i-l/2,j-\/2 


\(^~x^\h2/ ho)i,j (l~ x hih2 j fJ-o)i— l^j 

“I - (T~xh\ h2 / [J-o)i.,j — 1 (l~xhl h2/(J-o)i— 1 ,j — 1 ]/ (2 Ax) 


and 
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— \fr z h\h 2 / Ho)i,j (T z h\h2/ — l 


9(T z h\h 2 / Ho) 

dz i— 1/2, j — 1/2 

+ (r z hih 2 / Ho)i-i,j ~ (T z hih 2 /no)i-i,j-i]/(2&z) 

where Ax = (x z j — x*_ ij) and A z = (z Zi j — z z ,j- 1 ). Substituting into equation 4.4 and 
gathering terms: 


+ At(hi h 2 / Ho)i,j ( )i,j — 

+ (hi — ^2)*-ij-i + (hi — h 2 )i-\,j 4- (hi — h2)i,j-i] 

"k [(r x hih 2 / (i~xh\ h 2 J Ho)%— i,j (tj/ii/ij/ / io)jj_i](Af/ Ax) 

+ [( T zhih2/ Ho)i-l,j-\ + (T z hih 2 / ~ (T z hih 2 / Ho)i-l,j](At/ Az) 

Now solving for h2,»,j» the oil film thickness at the new time level: 

h2,i,j = {[^i,ij + (hi — h2)i~i,j-i + (hi — h 2 )i~i,j 4- (hi — h2)i.j-i] 

d - \(xxhih2 / Ho)i — i,j — i ~k (r^hi/^/ Ho)i — i,j (x x hih 2 / Ho)i,j — i](Af/ Ax) 

+ [( r zhih 2 / Ho)i-l,j-l + (T z hih 2 / Ho)i,j-\ — ( T zhih2/ Ho)i-\ j](A^/A2)} 

/[l + A t(h\/ Ho)i,j( ~^ + 

Equation 4.6 is the 2D Box-Implicit algebraic equation which allows us to numerically 
solve the 2D thin oil film equation, equation 4.1, at an interior node. The Box-Implicit 
node equation, as derived here, is second-order accurate in space, with “quasi” -higher-order 
treatment of the time variation. The equation is conservative of the oil mass. The Box 
method is known to be unconditionally stable. The related boundary condition treatment 
required to solve a problem is dealt with in the next subsection. 

The interior node equation, equation 4.6, clearly can be solved a single node at a time, 
with each interior node being solved sequentially. 

The form of the Box-Implicit interior node equation, as given by equation 4.6, is 
incomplete in that it does not include additional effects, for example, pressure gradient 
effects. If the oil film is sufficiently thin, these additional effects are of minimal significance. 
However, a more complete form of the interior node equation is: 


(4.5) 


(4.6) 


h%l,j ~{[hi,i,j + (hi — h2)i-\,j-\ 4- (hi — h 2 )i-i,j 4- (hi - h 2 )i,j~i] 

+ [( 7 ® — n x )j_i ; .,_i -t- (t x — n x )i_ij — (t x — n I ) i , J _i](A</Ax) 

~k [(Tz — n z )j_i j_i 4 - (T z — — (T z — n z )i_i j](Ai/ A z)} 


/[I 4 - At(— 


nl fc] 


Ax 


+ 


nL fc] 


Az 


)»,j] 


(4.7) 
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where 


T x =r x hih 2 /Ho = Tcos('y)hih 2 /Ho 
T z =T z h\h 2 / Ho = Tsin(i)hih 2 / Ho 
n z =(dP/dx - p 0 g x )h\h 2 {h\ 4- h 2 )/3p 0 
U z =(dP/dz - po9z)h\h 2 (h\ + /i 2 )/3p 0 
T x =r x h\/ Ho = T CQ${l)h\f Ho 
T z =r z h\/ Ho = rsin^hi/Ho 
ft ? 1 =(dP/dx - Po 9 x )h x {h x + h [ f~ l] )/3Ho 

nf ] =(dP/dz - + 4 fc_1] )/3/Xo 

JO] 

A: =1, fcmax = node iteration level 

To derive equation 4.7, we have made use of the numerical quadrature given by: 


f h?dt = ( — + — )hih 2 (t 2 - ti ) (4.9) 

Jti 2 

As with equation 4.3, equation 4.9 is first-order accurate in time, but is a quasi-higher- 
order time variation treatment since the numerical quadrature is exact for the practical 
case of h « t _1 . 

Due to the presence of h 2 in the II z and II 2 terms, iteration at each node (k = 1, fcmax) 
is required to advance the solution to the new time level. Also, the {i,j, etc.) subscripts 
are absorbed inside the parentheses as required to allow evaluation of each term of equation 
4.7 at the proper locations. 

An important property of the interior node equation (either equation 4.6 or 4.7) 
is that the ID self-similar problem solution for constant r (eq. 3.8, hrt = HoX) exactly 
satisfies either of these algebraic relations. This property of consistency of the interior node 
algebraic equations with the ID self-similar problem considerably enhances the accuracy 
of the solution method for most practical problems encountered. 
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Surface Tension Terms 


Equation 4.7, as shown, does not include surface tension effects, and, thus, the hyper- 
bolic nature of the partial differential equation is exploited by requiring only a node-by- 
node solution procedure. The inclusion of surface tension incorporates an elliptic feature 
and the solution molecule at node (i,j) should then include, for the additional surface 
tension terms, information from the nodes between ( i - 2 ,j) and ( i + 2 ,j) and between 
(*, j ~ 2) and (i,j -I- 2). 

The simplest surface tension thin oil film approach is to recognize that the time scale 
which characterizes surface tension is typically much greater than the time scale which 
characterizes a thin oil film acting under shear stress. Thus, the surface tension adjustment 
to the wall pressure at each node will be essentially constant during the integration time 
step, At. A sufficient treatment of surface tension for most thin oil film calculations would 
be an explicit approach where the adjustment to pressure is made at the known time level, 
ti, prior to each time step: 

P = Pair + Pa, where P a = -<r[(hi) xx + (h 1 ) zz ] (4.10) 

An example of an implicit elliptical surface tension treatment is given in the applica- 
tion section for problems where surface tension acts only in the x direction. However, we 
normally would take the above rather simplistic explicit approach due to the predominance 
of the viscous terms over surface tension for the thin films which we likely will have an 
interest in solving. This is done in order to maintain the simple to code node-by-node 
solution procedure. 
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Numerical Boundary Conditions 


The solution at the interior nodes must be started with an initial condition, where an 
initial value of oil film thickness h(x , z, to) is assigned, and then the solution is advanced to 
the next time level starting at boundary nodes. Thus, numerical treatments of the initial 
and boundary conditions are required. The types of boundary conditions required include 
a leading edge case, a corner leading edge case, a plane of symmetry case, a general surface 
streamline case and a boundary wall (no-flow) case. An important consideration for 2D 
oil film problems is that the boundary condition treatment will depend on whether the 
characteristics point into or out of the oil film domain at any particular location on the 
boundary. 

Leading Edge Boundary Condition 

The leading edge of an oil film is the contact line separating the region of the test 
surface covered by the oil film from the region of the test surface exposed directly to the 
aerodynamic flow. The wall shear stress at the leading edge will cause oil to flow from the 
leading edge into the region covered by the oil. A leading edge may be used as a boundary 
condition since the oil film height at the leading edge is known (h = 0) , and the direction of 
the wall shear stress, r x and r 2 , suggests that the convective velocities, U c and W c . within 
the oil film immediately adjacent to the leading edge allow solution at the interior nodes. 
In contrast, at a trailing edge, the wall shear stress will have a direction pointing out of 
the oil film region. The advancement of the oil film trailing edge over a fresh surface can 
be a difficult subject, involving finger instabilities, and is not treated in this present work. 
Therefore, the trailing edge is not a suitable boundary condition for a thin oil film solver. 

Over an extended period of time for which air flow occurs past a thin oil film, the 
leading edge of the thin oil film will eventually move in the direction of the aerodynamic 
flow, uncovering the test surface. However, here we consider the use of oil films in aerody- 
namics testing, and for moderate run times the oil film leading edge may be considered to 
be stationary. 

Further, a complete consideration of the fluid mechanics occurring at the leading 
edge of a thin oil film will include surface tension effects. At the leading edge, a contact 
line occurs which is defined by the juncture of the solid-air, solid-oil, oil-air interfaces as 
depicted in figure 5. Discussion of surface tension effects at the leading edge will be deferred 
to the applications section (Section 5). However, the amount of oil mass present in the 
control volume associated with the boundary nodes is quite small and for many practical 
cases, the simplified treatment of the leading edge presented here proves sufficient. 

The simplest leading edge treatment is to align the grid so that either the i = 0 or 
j — 0 boundary nodes line up precisely along the actual oil film leading edge. The oil film 
thickness for these boundary nodes is then simply held at h = 0, and the solution actually 
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starts at the next row of interior nodes, i = 1 or j = 1 , using a variant of the interior node 
equation 4.6: 


* 2,1 ,j 


— [*i,i,j + (* i - * 2 ) 1 , j-i + (h x h 2 At/p 0 ) ij-i(g^ ~ 
/[I + (hi At / + |^)iJ']’ j = 2,-,imax 


(4.11) 


Equation 4.11 is the boundary condition equation for the i = 1 row of interior node 
points for a leading edge at i = 0. Derivation of equation 4.11 is based on a control 
volume analysis which makes use of the mass flux in the z-direction integrated over the 
face between (l,j) and (0,j): 


rt 2 /*Ax 

/ / T z h\h 2 / ix 0 dxdt = (T z hih 2 / n 0 )i,jAtAx/2> (4.12) 

Jti Jo 

Note that at the leading edge the oil film thickness varies as h ^ x/t between nodes 
at * = 0 and t = 1. A similar node equation is easily derived for the case of a leading edge 
at j = 0. 

The intersection of two leading edges, with the i = 0 row of boundary node points 
forming one leading edge and the j ' = 0 row of boundary node points forming the other 
leading edge, creates the natural place, (i = 1, j = 1), to start the process of solving for 
interior nodes. Again, a control volume analysis for this corner node with /i 0 ,i — /ii,o = 
* 0,0 = 0, and assuming a linear variation of oil film thickness toward the node at (1, 1) 
leads to the corner leading edge node equation: 


* 2 , 1,1 — *i,i,i/[l + At(h\j Ho)i A (-^~ + ^)i,i] 


(4.13) 


Plane of Symmetry Boundary Condition 

A plane of symmetry for the aerodynamic flow may serve as a boundary condition 
for the thin oil film equation. Consideration of the plane of symmetry (aligned with the 
rr-axis) leads to the following relations: 

a. dh/dz = 0, but h ^ 0. 

b. 7 = 0. but dj jdz ^ 0, which implies r 2 = 0, but dr z /dz ^ 0. 

C- po9z = 0, but dp a g z jdz ^ 0. 

d. dP jdz — 0, but P ^ 0, and d 2 P/dz 2 ^ 0. 

The surface streamline angle, tan( 7 ) = t z /t x . may be derived from the known shear field. 

We assume the row of boundary nodes, (i = l,imax, and j = 0), is aligned along 
the line of symmetry formed in the thin oil film. The numerical form of the symmetry 
boundary conditions is based on equation A. 3, repeated here as: 
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dh drh 2 / 2p 0 

dt dx 


+ (rh 2 /2fj, 0 )-^ — 0 


Evaluating at (i — 1/2,0) using the same numerical techniques leading to the interior 
node equation 4.6 results in: 


1 1 dy 

^2,1,0 =[^l,z,0 + {h\ ~ ^ 2 ) 1 - 1 , 0 + At(T x h\h2 / ~~ 2 


1 1 dy 

/[I + At(r x hi / 2 dz 


(4.14) 


Equation 4.14 applies in the absence of pressure and/or gravitational effects. For the 
more general case where pressure and/or gravitational effects are included, we start with 
the complete form of the plane of symmetry equation: 


dh d ,r x h 2 r x h 2 d'y 

dt dx v 2p a ’ V 2 p 0 ’dz 


d u dP_ ^ 3 i ( h \( d2p 

dx [{ dx Po9x) 3p Q ] 3/x c dz 2 


dpo9z 

dz 


) = 0 (4.15) 


An equation similar to the interior node equation 4.7, but for the boundary nodes 
located on the plane of symmetry can now be derived from equation 4.15: 


^2,1,0 — [hl,i,0 + (hi — /l 2 )i-l ,0 + Af( - ; — 2^ X dz 

/[1 + A(( ^ + I fl |2 + fiW,,,] 

where 

T x =T x hih 2 /Ho = rcos^hxh-i/no 
n x =(dP/dx - p 0 <?x)[hih 2 (hi 4- h 2 )/3p 0 ] 
n n =(d 2 P/dz 2 - dp 0 g z /dz)[h-ih 2 {hi + h 2 )/6/z 0 ] 
T x =T x hi/p 0 - rcos( 7 )hi//i 0 
ft? 1 =(dP/dx - p 0 Qx)[hi(h\ + h [ t l] )l3po\ 
n[f ] =(d 2 P/dz 2 - dp 0 g z /dz)[hi(hi 4 h£~^)/6n 0 ] 

U [0] _ u 

k =1, fcmax = node iteration level 


(4.16) 


(4.17) 
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The dp 0 g z /dz term may appear if the plane of symmetry of the body cuts vertically 
through the body. The addition of the pressure and gravity terms requires an iteration 
at each node as was done for the interior node equation 4.7. Further, as with the interior 
node equation 4.7, if surface tension terms are to be included, we adjust the pressure using 
equation 4.10 in a global iteration scheme. 

Wall or No-Flow Boundary Condition 

A wall boundary condition may be used if the thin oil film is bounded on one side by 
a solid wall. The development is similar to the plane of symmetry boundary condition in 
that there is no oil flow through the boundary. Indeed, in the absence of pressure gradient 
or gravity effects, the plane of symmetry equations, equations A.3 and 4.14, may be used 
for a solid wall boundary. 

If the effects of pressure gradients, gravity, and/or surface tension are to be included, 
however, the dP/dz and dh/dz gradients need not be zero for the wall boundary conditions. 
As a consequence, a preferable form of the applicable governing differential equation for 
the thin oil film along a wall boundary is: 


dh d ( r s h 2 ( rj^_ .dy <L\r dP \ h3 i 

dt + dx { 2p 0 ) + [ 2p 0 ] dz dx [[ dx Po9x) 3 fxj 


h 3 , ( d 2 P dp o9z dP h?_dh 

3 n 0 )K dz* dz } { dz Po9z) p 0 dz 


= 0 


(4.18) 


The following numerical quadrature should prove useful in deriving the wall boundary 
node equation: 



, dh 


dhj 


h —dt - -h t fc 2 (^ + 


dh2 

dz 


)(*2 - * l ) 


(4.19) 


Since the wall boundary condition has not been tested in this present study, we defer 
in presenting the form of the wall node equation which includes the pressure gradient 
effects. 

Other potentially useful boundary conditions include a plane of symmetry intersecting 
a solid wall, a solid wall intersecting another solid wall, and a general surface streamline. 
The general surface streamline boundary condition can be derived for the situation where 
the oil streamline is known and the boundary is chosen so as to align with the known 
surface streamline. The plane of symmetry boundary condition is a special case of the 
general surface streamline boundary condition. The general surface streamline boundary 
condition should prove most useful for a mapped grid and is based on equation A.3. 
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Finite- Volume Upwind-Implicit Direct Solver: Interior Node 

An alternate approach to the finite-difference Box-Implicit numerical method is the 
Finite- Volume Upwind-Implicit numerical method. Here, we briefly develop only the in- 
terior node equation for a ID thin oil film case, with U c > 0. The extension to 2D and 
a change in sign of U c is not considered here. Further, the pressure gradient, gravity, 
centrifugal and surface tension terms are not included, since these terms may be treated 
in a manner similar to that employed for the Box-Implicit numerical method. 

Figure 6 shows a ID interior node computational molecule for the Finite- Volume 
Upwind-Implicit method. For the Finite- Volume method, hi represents the “average” oil 
film thickness within the node volume which lies between x* and Xi+i. When needed, we 
may use the approximation that the height hi is located at the midpoint of the node, (xi + 
Xi+i)/2. In contrast, for the finite-difference method, hi represents the oil film thickness 
at Xi and the average oil film thickness between Xj_ i and x* is given by (hi + hi- 1)/2. 

A mass balance during the time interval between t\ and t 2 for the control volume 
defined between Xi and Xi + i is: 

pt,2 f*t2 

( Am cv )i + ( / F i+ i /2 dt - / Fi_ 1/2 dt ) = 0 
Jt 1 Jti 

where / t * 2 F t dt = (p 0 r x h- l h 2 At/2fj, 0 ) i . Observe that the mass flux, F t . is evaluated at the 
cell volume midpoint, Xj +1 / 2 = (x^ + Xi+ 1)/2. Thus, the value for (r x )i is actually the 
value for the wall shear stress at x i+ i/ 2 . 

We do not directly know the mass flux at the control faces (i.e., F{ + 1/2)1 but rather 
we know the state of the fluid at the midpoint of the control volume (i.e., hi and F,). 
The mass flux at the control volume face is obtained by a second-order upwind biased 
extrapolation. 


Fi+1/2 = (3Fj — Fi_i)/2 (4.20) 

A first-order flux limit relation is substituted to eliminate overshoots and oscillations 
where dF/dx changes sign: 

Fi+x/2 = Fj, if (Fi - F,_ 1 )(F t _ 1 - F t _ 2 ) < 0 (4.21) 

For most nodes the flux limiter is not applied and the mass balance for node i becomes, 
with rearrangement: 


/•l L \ A _ , A£ [0/ r x hih 2 x A ( T xh\h 2 ^ ( T x h \Jv2 x ]_ n 
(h 2 — h\)iAx H — — 13( )i — 4( )i_i + ( )t— 2j — 0 

4 Po ho ho 


(4.22) 
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Rearranging, we obtain the ID Finite- Volume Upwind-Implicit interior node equation 
(without Flux Limiter): 


{h 2 )i = [(hi)i 


/ rj/ii/12 , 

Ax [i 0 1-1 4 Ax 


At t T x h\ h 2 ^ i 3 At ( T x h\ ^ 
V — )i-2\/[i + TTT:! - ~ id 


fj-o 


4Ax 


Vo 


(4.23) 


For comparison purposes, the ID Box-Implicit interior node equation (eq. 4.6) is: 


(M. = [(Mi + <*, - ft,).-. + £(^).-.]/[i + 


Numerical boundary conditions are developed similar to the Box-Implicit method, but 
the boundaries occur at the node faces, rather than where the node value for h is known. 
For example, for a leading edge at x = x 0 , h i=0 ^ 0 since (xi - x 0 )h i=0 represents the 
oil volume contained in the control volume associated with the i = 0 node. Rather, the 
leading edge condition implies that F_i/ 2 = 0 for the control volume face located at xo- 
Since, in the absence of surface tension, the partial differential equation is hyperbolic, 
the ith node depends on only those nodes upwind (e.g., i — 1, etc., for U c > 0) and the 
nodes may be solved sequentially. Also, notice that the Finite- Volume Upwind-Implicit 
method is mass conservative due to the consistent mass flux treatment at each face. The 
method is second-order accurate in space, with a quasi-higher-order time variation treat- 
ment. Further, the node equation 4.23 is consistent with the ID time self-similar solution 
for the case of constant r x . 

A favorable comparison of results from the Finite- Volume Upwind-Implicit method 
with results from the finite-difference Box-Implicit method will be made in the application 
section (Section 5) for a ID test problem. 


4.2 INVERSE NUMERICAL SOLUTIONS 


Numerical techniques for the inverse solution of the 2D thin oil film partial differential 
equation provide the basis for the experimental determination of the wall shear stress 
distribution generated on a 2D surface by a 3D flow. Because the inverse method described 
here is closely related to and derived from the Box-Implicit direct method described above, 
a brief derivation is provided of the interior node equation and boundary conditions. In 
the subsections below, we describe both a Two-Time-Level Box-Implicit inverse method 
and a One-Time-Level Box-Implicit method. The Two-Time-Level method is useful for 
experiments where sufficient optical access is available to view the test surface during the 
experiment, while the One-Time-Level method is suited to those experiments where an 
image of the oil film is acquired after the experimental wind tunnel run. 
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Two-Time-Level Box-Implicit Inverse Method 


To apply the numerical technique described here, the oil film thickness at two dis- 
tinct times, h{x,z,t\) and h(x, z, £2)1 over the entire test surface region of interest is as- 
sumed to be known from, for example, optical measurements 1-4,8,9-11 . Also required are 
measurements 9-11 of the surface streamline direction, 7 Or, z), over the same region. The 
aerodynamic flow over the surface is assumed to be steady between these two times. W hen 
pressure gradient effects are important, it is necessary to know the wall pressures during 
the same time interval, but the inclusion of pressure gradient, gravity or surface tension 
effects does not alter the solution strategy. 

In the context of understanding the solution of the 2D thin oil film equation for the 
wall shear stress, "7*, equation 2.2 may be rewritten as: 


d r x h 2 + d r z h 2 +s _q 
dx 2/z 0 dz 2 Ho 


(4.24) 


Where S' is a source term absorbing all those terms of equation 2.2 which do not 
contain the wall shear stress, r. 

Next, apply a coordinate rotation as given in Appendix A (see eq. A. 3) to equation 
4.24 by the angle 7 from the (x, z) coordinate system to the (s, z) coordinate system aligned 
locally with the wall shear stress: 


d rh 2 rh 2 d'y ^ 

7TT~ + +S = 0 

os 2 Ho 2[z 0 oz 


(4.25) 


By examining equation 4.25, notice that we now have a first-order differential equation 
for t that can be solved along the s-direction which lies aligned with the characteristic 
direction, 7. The form of equation 4.25 emphasizes the hyperbolic nature of the thin 
oil film equation when solving for r, and that the characteristic direction for the inverse 
solver is given locally by 7, regardless of the nature of the source terms (including possible 
pressure gradient, gravity or surface tension effects). In contrast, when equation 2.2 is 
solved in the direct mode for h , inclusion of the pressure gradient and gravity terms can 
significantly affect the characteristic direction (U C ,W C ). Further, if surface tension terms 
are included when solving in the direct mode for h , the nature of the equation changes 
from hyperbolic to elliptic. The consistent hyperbolic (upwind) nature of the thin oil film 
equation when solving for r enables the inverse solver to incorporate the pressure gradient, 
gravity and surface tension effects in a more straightforward manner (without iteration) 
than the direct solver. 

Although a possible solution strategy would be to solve equation 4.25 along identifiable 
surface streamlines, the approach here is to solve equation 2.2 on a 2D grid in a manner 
related to the direct solver described earlier. Such an approach is easier to implement for 
the general experimental case. A particular advantage to the chosen approach is that the 
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difficulty alluded to in Appendix A of finding the s-n coordinate system is thereby avoided. 
Knowledge of 7, the surface streamline direction, allows us to define the proper domain of 
influence, and, thus, which grid points are required to be included in the solution molecule. 

The same solution molecule used to describe the interior node for the direct solver 
is also used to describe the inverse solver (see fig. 4). The domain of influence shown in 
figure 4 assumes 0<7<7r/2. A rectangular 2D grid is assumed, with even spacing of Ax 
and A 2. Such a 2D array of grid points might, for example, correspond to the 2D pixel 
arrangement for experimental data obtained from a series of digital camera images of a 
thin oil film on the test surface. 

The control volume analysis which leads to the Box-Implicit direct solver interior 
node equation still applies for the Box-Implicit inverse solver. Thus, the interior node 
equation which forms the basis of the 2D inverse thin oil film solver is a straightforward 
rearrangement of the interior node equation for the Box-Implicit direct solver, equation 
4.6: 


( r )i,j — {[(hi — ^2 )ij + (hi — h.2)t-l,_j-l + (h\ — h 
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(4.26) 


Although the size of the time step between the two images, (£2 — i 1 ) , may influence 
the accuracy of a solution, the stability of the inverse solver, which marches in space, is 
not dependent on the time step size. 

A more complete form of equation 4.26 may be derived by rearrangement of equation 
4.7 (rather than eq. 4.6) so as to include the pressure gradient, gravity and/or surface 
tension effects. The inverse solver, even with the inclusion of these terms, does not require 
either global or node iteration, unlike the direct solver. 

Boundary conditions for the inverse solver may be derived from related boundary con- 
ditions for the direct solver, generally by rearrangement. For a leading edge, the boundary 
condition for the inverse solver is analogous to equation 4.11: 


2 T T 

r i ,j = [(hi — h2)i,j -I- (hi — h2)ij_i + (hih2At/ fx 0 )ij-i(--^~ — •^)i, J _i 
/[(hih 2 At/fi 0 )ij(^^- + 3 = 2, ..,jmax 

For a plane of symmetry, the boundary condition for the inverse solver is analogous 
to equation 4.14: 


(4.27) 
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(4.28) 


1 1 d^f 

r i,o = [(h\ ~ ^2)1,0 + (^1 — ^2)1-1, 0 + {rhykiAt/ /x 0 )i-i,o( ^^ — ^ ^)<-i,o] 

1 1 dy 

/[(hih 2 At/fj, o )i' 0 (— + )i,o] 


The inverse solver then solves on a node-by-node basis, starting with the boundary 
nodes and then proceeding to each of the interior nodes as the information for each neigh- 
boring upwind node becomes complete. For example, provided 0 < 7 < 7r/2, the nodes 
would be swept according to the following simple row-column strategy: 

a. For £=1, and j = 0, advance the plane of symmetry /leading edge boundary node. 

b. For i = 1, and j = 1, jmax, advance each node in the leading edge boundary. 

c. For j = 0, and i = 2, imax, advance each node in the plane of symmetry boundary. 

d. For an outer loop of i = 2, imax, and an inner loop of j = 1, jmax, advance each of the 
interior nodes. 

The rectangular 2D grid array described above is a somewhat limiting feature of 
the present treatment of thin oil films. A more general grid mapping transformation of 
irregularly spaced nodes located at (xi,j, Zij ) to evenly spaced coordinate system (Qj, rjij) 
may be accomplished, but is not treated here. 


One-Time-Level Box-Implicit Inverse Method 


Optical access to the test surface may be limited during an actual wind tunnel run. 
Thus, some researchers prefer to acquire a single image of the oil film taken immediately 
at the finish of the wind tunnel run. In this subsection, we derive an inversion method 
suitable for analyzing such a One-Time-Level 2D oil film thickness distribution. 

The wind tunnel is assumed to have run sufficiently long and the oil film is thin enough 
that pressure gradient, gravity and surface tension effects are negligible. Further assume 
the oil film distribution has achieved time self-similarity, where h(x,z,t ) = H(x,z)/t. 
Thus, we start our derivation with equation 3.3, rewritten here as: 


8 ,t x H 2 
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Wind tunnels seldom start immediately at the desired run conditions. Presumably, 
the oil film is applied before the wind tunnel starts. Some initial time must occur to 
establish the desired run condition, the tunnel then is held at the desired run conditions 
for a period of time, and then the tunnel takes a small amount of time to shut down. 
During the entire time of the oil flow, the wall shear stress varies, as does the dynamic 
head, q x . Following the suggestion of Monson, Mateer, and Menter 7,8 , it seems best to 
consider that Cf = r/q^ is likely to vary less during the time of the oil flow than will r. 
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We, therefore, rewrite the 2D self-similar equation above as: 


+ £«?„*•) - 2 VjT' = 0 (4.29) 

where Cf = r/q nom , C/ x = C/C 0 S 7 , C/ z = C/sin 7 , and g nom is the nominal or desired 
wind tunnel dynamic head. The measured wind tunnel dynamic head during the wind 
tunnel run varies and is given by <joo(t). Additionally, the oil viscosity is temperature 
sensitive and may vary during the wind tunnel run, which requires that the oil viscosity 
variation with wind tunnel time be determined, fj, 0 (t). This can be done by measuring the 
test surface temperature and then using a temperature calibration for the oil. 

The derivation of equation 4.29 from equation 3.3 is not strictly rigorous in that we 
integrate equation 3.3 over the time interval from t s to if, and approximate the integral: 

[ qoo(t)H 2 /fj, 0 (t)dt « H 2 [ qoo(t)/n 0 {t)dt 

Jt s Jt s 

The difficulty with this is that the oil will not initially have been time self-similar, and 
H = h/t will actually vary with time. For the purposes of this single time level scheme, 
however, we consider this error source acceptable. The success of the overall method will 
be considered in the applications section. 

A One-Time-Level Box-Implicit numerical inverse equation form of equation 4.29 suit- 
able for interior nodes may be derived in a manner similar to equation 4.26 which solves 
the interior node for the Two-Time-Level numerical inverse scheme. The One-Time-Level 
interior node equation is: 
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(4.30) 


where 0=1/ q^t)/ fj, 0 (t)dt. Further, because the oil film distribution for only one 
time level is used, the time level subscript on h is dropped. 

Leading edge and plane of symmetry numerical boundary conditions for the One- 
Time-Level Box-Implicit method may be derived in a manner similar to that of equations 
4.27 and 4.28. 
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5. APPLICATIONS 


An understanding of the behavior of thin viscous oil films can be achieved through 
the study of selected ID and 2D model problems. Selection of these model problems 
is made so as to emphasize those physical, mathematical and numerical features which 
might arise in considering the use of the methods described in this paper. The relative 
influence of wall shear stress, wall pressure, gravity, centrifugal and surface tension effects 
can be determined. Comparison between numerical solutions, both direct and inverse, for 
those ID and 2D thin oil films for which closed form solutions can be found will reveal 
the suitability of the numerical methods for more general problems. The accuracy of 
approximate solution methods used by others in the analysis of experimental data can 
be assessed. Demonstration of the rigor, accuracy, and ease-of-use of the mathematical 
treatment of these thin oil films, as described here, tends to reinforce confidence in the 
continued development and expanding use of thin oil films in the important measurement 
of wall shear stress. The utility of the direct and inverse solvers for practical applications 
should become apparent from results of those cases presented. 

Variable Wall Shear, ID Case 

The first problem considered is a simple ID model problem where wall shear stress 
varies linearly. For this problem, we assume the pressure gradient is either zero or negligi- 
ble. Gravity and surface tension effects are also ignored. The solution is known from the 
ID analytical relations and we may evaluate the accuracy of both the Box-Implicit direct 
and inverse solvers for this case, as well as alternate solution methods which have been 
used in other studies for analyzing experimental data. 

Assume the wall shear stress varies as r(N/m 2 ) = a + bx = 20 + lOOx, where x is 
in meters, with the leading edge of the thin oil film at x = 0 and the film extending to 
x = 0.1 meter. Assume the oil film has a kinematic viscosity of v 0 = 100 centiStokes, with 
a density of p a = 1000 Kg/m 3 . These properties are similar but not identical to commonly 
used silicone oil. For the numerical solutions, assume an initial oil film thickness at t = 0 
of ho = 10 microns and a time step of dt = 0.1 second. 

A closed form time self-similar solution can be obtained for this case using equation 
3.7, the analytical relation for a ID thin oil film: 

H — ht — (po/r) l/2 jjpo/rY^dx = 2^(1 - (1 ~ bx/a) i/2 ) ^ 

Figure 7a shows the oil film thickness at several times (t = 0, 20, 40, 60. 80 and 
100 seconds) from the analytical solution above, and from the ID version of the Box- 
Implicit direct solver. An evenly spaced grid of 100 points is used for the numerical solver. 
The analytical solution assumes an infinite initial film thickness. Some of the observed 
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differences between the analytical and numerical solutions are due to this differing initial 
condition. However, very good agreement (< 1%) between the analytical and numerical 
solutions can be seen to occur for those regions and times where the oil film has reduced 
in thickness to 3/4 or less than the initial thickness. We consider that the decrease in 
differences between the numerical solution and the analytical solution appear to be more 
closely associated with the decrease in the oil film thickness relative to the initial thickness 
rather than with the increase in time. Clearly, however, the numerical solution approaches 
the analytical time self-similar solution as the oil film thins. 

Another observation from figure 7a is that for “small” times a corner exists in the 
oil film thickness distribution (e.g., x 0.04m for t = 20 seconds). For larger times, 
the corner has convected out of the solution domain. Note that by integration, for the 
dh/dx = 0 region with r = a + bx, of the ID form of equation 3.1 for the time variation of 
h , we obtain: 


h — ho/(l + hobt/2/j, ) (5.2) 

A piecewise analytical solution can thus be formed for this problem with the smaller h 
from equation 5.1 or 5.2 being selected. The observed corner occurs where equations 5.1 
and 5.2 intersect, x c = (a/6)[(l + hobt/2fj,) 2 - 1]. 

For the Box-Implicit direct solver, an oscillation in the solution appears to originate 
from this corner. This oscillation is numerical rather than physical. To solve this numerical 
oscillation, a “Flux-Limit” concept is adapted to the Box-Implicit algorithm. For this 
Flux-Limit concept to be applied, a first-order (rather than second-order) form of the 
Box-Implicit algorithm is used for those nodes where the slope of the oil film thickness 
changes sign. The solution from the Box-Implicit Flux-Limit algorithm is shown in figure 
7b. Examination of the solution in figure 7b at, for example, t = 20 seconds shows that the 
oil film slope changes sign at only two locations, meaning the first-order treatment at these 
two nodes is sufficient to eliminate the corner oscillation. Clearly, the corner oscillation is 
eliminated with excellent agreement otherwise. 

Figure 7c shows the solution from the Finite- Volume Upwind-Implicit direct solver 
algorithm. Although the solution appears to be somewhat better behaved than the Box- 
Implicit direct solver, a corner oscillation still appears. A Flux-Limit form of the Finite- 
Volume Upwind-Implicit algorithm eliminates this oscillation as is evident in the solution 
given in figure 7d. 

An adaptation of the MacCormack 22 explicit algorithm to the thin oil film problem 
leads to the solution given in figure 7e. The explicit algorithm is stable since the effective 
stability index (CFL= U c At/ Ax ~ rhAt/2(i 0 Ax) for this problem is less than 1 . No corner 
oscillations occur, although slight differences (< 2%)in the oil film thickness distribution 
for t = 20 and t = 40 seconds compared to the other numerical solutions are observed. 

A further test of the numerical methods is to consider the accuracy of the inverse 
solver. To accomplish this, we use the analytical oil film thickness for two times, t — 80 
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and t = 100 seconds, as the known input to the inverse or r numerical solver. The 
r vs x solution from the inverse solver can then be compared with the original known 
r = 20 4- lOOx distribution. Figure 8a shows a r vs x distribution from the Box-Implicit 
inverse solver. The comparison with the known result is excellent, within 1%. A further 
test is to use the oil film thickness distribution from the Box-Implicit direct solver as the 
input to the Box-Implicit inverse solver as shown in figure 8b. Clearly demonstrated is the 
accuracy of these methods for both direct and inverse solutions for this test case. 

Shown in figure 8c is the r vs x distributions derived by means of the One-Time-Level 
Box-Implicit inverse algorithm from the Box-Implicit Flux-Limit direct solver’s oil film 
thickness at times t = 40 and t = 100 seconds. Agreement within 1% occurs except for 
where the corner region of the t = 40 second Box-Implicit oil film thickness distribution has 
not yet convected out of the region of interest. The One-Time-Level Box-Implicit inverse 
algorithm requires the input oil film distribution to have become fully time self-similar 
(H(x) = /i(x, t)t ^constant), which has not occurred as yet for the t = 40 second numerical 
oil film distributions. A reasonable estimate for the time at which an oil film becomes fully 
time self-similar is given by t s ~ L/U c ~ Lujrho. For this problem, L = 0.1 meter, 
ho — 10 microns and t s ~ 50 seconds. When the input oil film thickness distribution has 
evolved to the point where it is time self-similar, the One-Time-Level Box-Implicit inverse 
algorithm yields quite accurate results. 

An alternate method of solving for r vs x that appears in the oil film literature is 
based on the ad hoc assumption that the local wall shear stress is inversely related to the 
local slope of the oil film height (only one time level is used): 


T5s( t )/( f ) (53) 

Figure 8d presents the result of application of this local slope method to determine the 
r vs x distribution from the analytical oil film thickness at t = 100 seconds. Comparison 
with the known r = 20 + lOOx distribution reveals obviously large errors (> 50% of the 
rise in r) which result from the use of equation 5.3. These errors can be shown to arise 
where dr/dx ^ 0. In particular, by rearrangement of equation 3.5: 


T 


[<t 


(5.4) 


The error in equation 5.3 is given by the (ht/2n 0 )dr/dx term. The increase in r of 
lOAT/m 2 over the region of the oil flow from x = 0 to x = 0.1m leads to an error of 6Af/m 2 
or 60% for this algorithm. It is strongly recommended that the relation of equation 5.3 not 
be used except for regions very close to the leading edge of the oil film (as was correctly 
done by Monson, Mateer, and Menter 7 ’ 8 ). Where the oil film is time self-similar and a 
single image is available, the One-Time-Level Box-Implicit inverse algorithm described 
earlier is recommended instead due to its ease-of-use and high accuracy. 
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Constant Wall Shear, Increasing then Decreasing Film Thickness, ID Case 

A feature related to the nonlinearity of the thin oil film equation is explored with 
the next problem. The study of the behavior and solution of nonlinear partial differential 
equations (PDEs) invariably leads to Burgers’ equation. A point to be made is that, where 
r and p Q are held constant, the simplest form of the ID thin oil film equation is identical 
to the inviscid form of the Burgers equation. Much of what has been learned from the 
study of Burgers’ equation about the behavior and solution of nonlinear PDEs directly 
applies to the present thin oil film work. In turn, we realize that the study of the thin oil 
film provides a quite practical and possibly interesting physical manifestation of Burgers’ 
equation. A certain irony exists in the realization that each oil flow study performed by 
an experimentalist is also actually a Burgers’ equation experiment. 

One CFD problem area studied using Burgers’ equation is the inviscid ID (or normal) 
shock. The present problem examines this shock-like behavior that can arise from the 
nonlinear term, even for a thin oil film. 

The ID form of the thin oil film equation is given by: 

dh drh 2 /2p 0 

dt dx 

An advantage of the thin oil film form is that we may form model problems where r 
varies and even changes sign, allowing us to examine a greater range of fluid physics. Such 
a model problem is examined in the next subsection. For the particular model problem 
under consideration, however, we assume that r = 25 N/m 2 , p = 1000Kg/m 3 , v = 100 
centiStokes (e.g., i v = 100- 10 -6 m 2 /sec) and p = pu = O.lKg/m-sec and that these remain 
constant. 

First note that oil thickness, h, is being convected in the mean at a velocity, U c = 
rhj 2p, which leads to the nonlinear h 2 term in the PDE. The characteristic wave velocity, 
found by expanding the d/dx term, is given by U w = rh/p, which is also equal to the oil 
velocity at the air /oil film interface at y — h. 

Further assume for this current problem that the initial (to = 0) oil film thickness 
increases linearly with r, followed by region of constant thickness, then decreases linearly 
with x followed by another region of constant thickness according to: 

h(t 0 = 0) =10 _5 x/ 0.04, for x < x a = 0.04m 

h(to = 0) =/io = 10 _5 m, for x a = 0.04m < x < Xb = 0.05m 

h(to = 0) =10 -5 (3.5 — x/0.02), for Xb - 0.05m < x < x c = 0.068m 

h(to — 0) =h\ = 10~ 6 m, for x > x c = 0.068m 

Figure 9a shows this initial (to = 0) oil film thickness distribution, as well as subse- 
quent development of the oil film for times of t = 4 , 8, 12, 16 and 20 seconds. The corners 
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initially at x a = 0.04m, Xb = 0.05m and x c = 0.068m are also identified in figure 9a to aid 
in describing the subsequent oil film development. 

By means of characteristic and mass-conservation based arguments, the development 
of the oil film thickness distribution with time can be established. Such piecewise analytical 
solutions of the oil film provide a basis for comparison with numerical solutions, as well as 
provide revealing insight into the shock-like behavior exhibited. 

Consider the region between xq — 0 and x a , where the oil film thickness increases 
linearly with x. Note that the convective velocity also increases linearly in x and the film 
tends to spread, leading to a decrease in the slope of the oil film thickness. Because the 
region between x = 0 and x a is linear with constant r and n Q , we know the ID time 
self-similar solution applies, where rh(t - 16) = /jl 0 x. Rearranging slightly, we find the 
slope for this region is Ah/ Ax = n 0 /x(t — 16). The corner at x a , where the linear region 
intersects the h = ho plateau, thus moves to the right at a constant wave velocity such 
that: 


x a (t) ~ x a (to) = ho/ (Ah/ Ax) = (rh 0 /fJ- o )(t - t 0 ) 

The wave velocity of the corner at x a is rh 0 / (x Q = 2.5mm/sec. Similarly, the corner 
at Xb is an artifact that moves at the constant wave velocity, rhoj ii 0 — 2.5mm/sec. Also, 
the corner at x c moves at the constant wave velocity, rh\/n 0 — 0.25mm/sec. 

The plateau between x a and Xb will retain constant oil film thickness since dh/dt = 
— (rhj fj,)dh/dx = 0, with corners at x a (t) and Xb(t) which appear to translate to the right 
at a constant wave velocity of rho/n- 

In contrast, between Xb and x c . the oil film thickness decreases linearly with x. Thus, 
the wave velocity, rh//x 0 , also decreases in x. The corner at Xb tends to catch up with 
the slower moving corner at x c and the oil film thickness slope becomes steeper, until 
a discrete jump of h 0 - hi occurs with Xb(t) = x c (t). The discrete jump forms when 
Xb(t ) = x c (t ) = 0.07m at t=8 seconds. 

Piecewise solutions for the oil film thickness at t = 4 and t = 8 seconds based on the 
above arguments are given in figure 9a. 

For times past t = 8 seconds, the discrete jump moves with a wave velocity associated 
with the average height of the jump, r(ho + h\)/2fi 0 = 1.375mm/sec. This wave velocity 
can be determined by considering a small control volume of length, Ax, with the jump 
just entering. The jump is from a plateau of constant thickness, ho, to another plateau 
of constant thickness. hi. Thus, the flow of oil into the control volume will occur at a 
constant volumetric rate, rho/2/x 0 . The flow of oil out of the control volume also occurs at 
a constant volumetric rate, tY/\/2\i 0 . The net rate of increase of oil volume in the control 
volume will be t(/iq — h\)/ 2(i 0 . However, the initial oil volume present as the jump just 
enters the control volume is hi Ax. The oil volume present as the jump just exits the 
control volume is ho Ax. The net change of the oil volume in the control volume as the 
jump travels the distance Ax is (ho - hi) Ax. The amount of time, then, for the jump to 
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travel through the control volume is At = ((/i 0 -/ii)Ax)/(r(ho-/if)/2^ 0 ). Rearrangement 
gives the wave speed of the jump: 

Uj = Ax/ At = r(ho 4- hx)/2n 0 

The piecewise solution for the time t = 12 seconds, as shown in figure 9a, results from 
application of the above arguments. 

When t = 15-^y seconds, the corner at x a will just catch up with the jump at x a — 
Xj = 78-^-mm. Subsequently, the ho plateau will no longer exist, and the average height 
of the jump will decrease as the jump moves to the right, and the wave velocity of the 
jump will also decrease. However, a consideration of the known time variation of the slope 
of the linearly increasing region (xo to Xj), along with the known time variation of the 
total oil volume, V = (0.431 — 0.0005f)10 _6 m 2 , in the domain (i = 0 to i = 0.1mm) 
allows calculation of the jump movement and oil film thickness distribution past t = 15fy 
seconds. The oil film thickness distributions for t = 16 and 20 seconds are also shown in 
figure 9a. 

Precise knowledge of the piecewise analytical solution for this “shock-like” ID model 
problem allows for assessment of the accuracy of the Box-Implicit method and the Finite- 
Volume Upwind-Implicit method for direct numerical solvers. These numerical solutions 
are given in figures 9b, 9c, 9d and 9e. 

With the exception of the region about the shock-like jump, the Box-Implicit and 
Finite- Volume Upwind-Implicit methods provide quite acceptable solutions. The jump 
region, however, leads to large oscillations in the numerical solutions particularly for the 
Box-Implicit method without the Flux-Limit treatment. For both the Box-Implicit and 
Finite-Volume Upwind Implicit methods, the Flux-Limit treatment improves the solution 
about the jump region, but with an apparent increase in jump wave velocity. This increase 
in jump wave velocity is likely because the rather simple Flux-Limit treatment is not 
mass-conservative as implemented. 

Another criticism of the numerical solutions obtained is that the ho plateau is not 
preserved, becoming rounded at the x a and Xb edges. Upwind methods are known to be 
dissipative and this appears to be the reason, with the Box-Implicit method being better 
behaved in this regard. 

2D Surface Tension Bubble Problem 

The next case considered is designed to demonstrate the inclusion of surface tension 
terms. The primary focus of the present paper is on the response of an oil film to a 
wall shear stress. However, in the vicinity of separation, the surface tension term can 
dominate. Thus, it is desirable to demonstrate that the surface tension terms are correctly 
implemented in the overall numerical procedure. To do this, we consider a problem where 
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the oil film is subjected to only the surface tension term, and, furthermore, the problem is 
designed so that features of the oil response are analytically known. 

In constructing the present problem, we do not solve the surface tension problem in 
general, but only demonstrate inclusion of the surface tension terms and also prepare for 
the 3D separation problem of the next subsection. We consider here, then, the simpler 
case where the oil film thickness varies only in the x direction. 

Where only surface tension effects occur, and then only in the x direction, equation 
2.2 becomes: 


dh d t h z .d, n 
dt + te { 3^ 0 [ fa {<Thxx)]} - 

Assume that oil is spread initially ( t = 0) so as to assume the shape: 
h(x , t = 0) = ho[l — (x/xo) 4 ], for — xq < x < xq 


(5.5) 


(5.6) 


This initial bubble shape is not in equilibrium and the initial rate at which the bubble 
changes shape is known by simple substitution of equation 5.6 into equation 5.5: 


Oh _ d h 3 d 
( dt )t=0 ~ dx { 3» 0 [ dx {ahxx)]} 

= 8-^^[l - (x/x 0 ) 4 ] 2 [l - 13(x/x 0 ) 4 ] 

VoX o 


(5.7) 


The bubble shape continues to change with time until steady state is reached. The 
steady state bubble shape will be given by: 


h(x, t = oo) = 1.2ho[l — (x/x 0 ) 2 ], for — xo < x < xq (5-8) 

Note the steady state bubble shape given by equation 5.8 satisfies equation 5.5 with 
dh/dt = 0, and also has the same oil volume as the initial bubble shape equation 5.6: 

/ Xo 

hdx = 1.6hoXo 

-XO 

An accurate, mass-conservative numerical direct solution method given the initial 
bubble shape of equation 5.6 should lead to the steady state bubble shape of equation 5.8 
and have an initial rate of change of oil film thickness given by equation 5.7. 

As can be seen from equation 5.5, the surface tension term is highly nonlinear ( h 3 ) 
involving a fourth-order derivative in space. To implement the surface tension term into 
the numerical solver, we choose here to use the linearized line-implicit method described 
below. Although the method is incorporated in the current example program for nodes 
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along the plane of symmetry as well as for purely interior nodes, we describe only the 
interior node treatment here. 

First, integrate equation 5.5 with time between times ti and £ 2 : 

h 2 ~ + ^[( 3 ^-) J t h z h xxx dt) = 0 (5.9) 

We linearize the h 3 h xxx term at £2 as: 

(h h xxx )2 ~ hj(/l2)xxx “h 3/l| /l2 (hi ) xxx 3/l| (hi )xxx 

Substituting into equation 5.9: 

h 2 — hi + ^[(-^— )(hi(/i 2 )xxx + 3hi/i2(hi) zxx — 2hl(hi) xxx )] = 0 (5.10) 

Now evaluate equation 5.10 at the midpoint, ( i — 1 / 2 ): 


h2,i — h\,i + /l2,i-i 


~ + ( 


a At 
3p 0 Ax 


XXX + 3 hih 2 (h\) 

XXX 2h?(hi) 

xxx]i 

— [hi (/i2) xxx + 3hi h2{h\)xxx 2hj (hi) ZXI ]j_i } = 0 


(5.11) 


where Ax = x l — x,_ 1 . 

The derivative, h xxx , may be found numerically by: 


( h xxx )i — Q{hi + 2 4- /3ihj + i + ejhj + 7ihi_i + 5jhj_2 

where the coefficients, (a, /?, 7 , 6 and e), are given in Appendix B. 

Equation 5.10, with substitution and rearrangement, takes on the following matrix 
form: 


A2ih.2,i-3 + Alj/i2,i — 2 + A0ih,2 y i- 


1 +. B0{h2,i + C0i/l2,i+l + Cljh2,i+2 — DOi 


(5.12) 


where 



C2i =0 

Ch =( 
COi =( 


aAt 3 


3fi 0 Ax' 
aAt 
3 fj, 0 Ax 




BOi — 1 + (— — )[*i,« e * ~ *i,i-iA-i 


AOi =1 4 ( 


3n 0 Ax' 

4 3/lii(a</li ) i+2 + A*l,i+1 + £4*1,4 + 7t*l,i-l 4 A*l,i-2)] 

vr>3 _ ^3 


3fi 0 Ax 


- *i,4-i e i-i 


— 3/lj 4- A— l*l,i 4 4- 7i-l*l,i-2 + <5i-l*l,i-3)] 


Ali =( 


aAt 

3fj. 0 Ax 


){h\A - Ai,i_i7i-i) 


Z)0i =/ii,i + 

+ 2h\ i {otih\^iJ r 2 + A *1,4+1 + £t*l,t + 7t*l,i-l + A*i,i— 2 ) 

*l,i+l 4 A- 1*1,4 + 4 7i-l*l,i-2 4 5i_i/li ? i_3) 


The matrix equation 5.12 is then solved by means of a scalar septa-diagonal solver 
written for the present work. The scalar septa-diagonal solver is quite similar to commonly 
used scalar tridiagonal solvers available, but is seven elements wide rather than only three. 

The initial oil film thickness distribution for the model surface tension problem, equa- 
tion 5.6 (with 100 centiStoke oil of 1.0 specific gravity and a = 21.010~ 3 N/m, ho = 10.0 
microns and xo = 0.1 meter) is shown in figure 10. The known analytical steady state 
oil film thickness distribution, equation 5.8, is also compared in figure 10 with the present 
steady state numerical results using the method described above. The excellent agree- 
ment of the numerical and analytical steady state oil film thickness distribution is clearly 
indicated. 

Figure 11 shows the initial rate of change in the oil film thickness, (dh/dt) t = o- for the 
model surface tension problem. Again the excellent agreement of the numerical results and 
the known analytical form (eq. 5.7) of this feature is clearly indicated. Suggested by this 
level of agreement for these features between the numerical result and the known analytical 
forms is that the surface tension terms are correctly implemented into the present oil film 
solver using the methods described in this subsection. 
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Saddle on Plane of Symmetry, 2D Case 


As a final example problem, the important practical case of a saddle located on a plane 
of symmetry is next examined. Such a flow occurs, for example, in the region ahead of a 
cylinder (or airfoil) mounted normal to a wall (or fuselage). Both a saddle of separation 
and a saddle of attachment are considered. These saddle cases are calculated by the 3D 
direct solver numerical methods described previously and are compared to the known closed 
form solution, also described previously. The additional effects of wall pressure gradients 
and surface tension for these saddle cases are also examined. This model problem is an 
example of a highly 3D separated flow, and the ability of the numerical methods described 
to accurately calculate this category of flows is a prime focus of the present study. 

For both saddle cases, we assume the saddle (where both r x and r z equal zero) is 
located at x = 0 and z = 0. The oil is applied initially ( t = 0) as a thin film of ho = 10 
micron thickness with a leading edge located at x = -0.01 meter. The oil properties were 
chosen to be nearly (but not identical to) those of silicone oil with a constant kinematic 
viscosity of v Q = 100 centiStokes, a density of p Q — 1000 Kg/m 3 , and a surface tension of 
a = 21 - 10~ 3 N/m. 

The wall shear stress is assumed to vary according to: 


r x = ax, and t z = bz 


For the saddle of separation case, a = -2000N/m 3 and b = 1500N/m 3 . 

For the saddle of attachment case, a — — 2000N/m 3 and b = 2500N/m 3 . 

In the absence of surface tension and pressure gradient effects, a time self-similar 
analytical solution is given by equation 3.21, repeated here: 


/ X (g + 6) 

h = Ht = (— — r)(l - ( — ) - 2a ), for 0 < x < x 0 
a + o xq 


(5.13) 


Further, the shape of the surface streamlines (passing through the arbitrary point 
(x s ,z s )) are given from the analytical solution by equation 3.19 repeated here: 


— = (— ) 6 /° = 
z s y x s y X 


(5.14) 


Figures 3a and 3b show the shape of the surface streamlines for the two saddle cases 
presently considered. 


Direct Solver, With and Without Surface Tension Terms 


Numerical solutions for the oil film thickness distribution were obtained using the Box- 
Implicit method for the two saddle cases, both with and without surface tension. Figure 12 
shows the numerical solutions at a time of t = 100 seconds, both with and without surface 
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tension. Also shown in figure 12 is the comparable analytical solution (eq. 5.13) for each 
case. Figure 13 is similar to figure 12, but shows the region near the saddle point in more 
detail. For this problem, the oil film thickness varies only with x and not z, even though 
the shear field varies in (x.z). However, the numerical solution procedures are indeed fully 
3D with h(x,z). Only the plane of symmetry solution is shown, but the z / 0 numerical 
solutions prove to be identical as would be expected for this particular test problem. 

The excellent agreement in figures 12 and 13 between the analytical solution and the 
numerical solution (without surface tension) gives a clear indication of the suitability of the 
Box-Implicit numerical methods we have implemented in the present direct solver. Only 
for the grid point located at x = 0 does a considerable difference between the analytical 
and numerical thickness solutions exist. The analytical solutions for both a saddle of 
separation and a saddle of attachment are cusped at x = 0. Furthermore, for a saddle of 
separation, the cusp at x = 0 leads to an infinite oil film thickness, h(x, 0) — » oo as x — t 0. 
The numerical solutions have some difficulty in resolving this cusp for the last grid point 
at x = 0. A further observation is that as the numerical grid becomes finer (not shown), 
the cusp in the numerical solutions at x = 0 becomes higher appearing to approach the 
analytical result. 

This sharp cusp, however, will not occur in a real oil film due to surface tension. 
Therefore, the inability of the numerics to resolve the sharp infinite cusp at x = 0 is of 
little practical significance. Also shown in figure 13 are the numerical solutions for t = 100 
seconds where surface tension is included using the line-implicit numerical procedures 
described and tested in the previous (surface-tension-only) problem. As can be seen in 
figure 13, the surface tension removes the sharp cusp which otherwise appears in the 
analytical and numerical solutions. Also, grid refinement (not shown) of the numerical 
solutions with surface tension no longer noticeably affects the oil film solutions obtained 
in the region of x = 0. 

Surface tension appears to more greatly affect the saddle of separation solutions than 
the saddle of attachment solutions, due to the more pronounced cusp occurring for the 
saddle of separation. From an order-of-magnitude analysis, the extent of domain affected 
by surface tension can be approximated by: 

= (5.15) 

M 

For the present saddle of separation test case, L a ~ 110 microns, while for the present 
saddle of attachment test case, L a ^ 80 microns. Based on figure 13, the observed domain 
actually affected by surface tension appears to be in the region of Arc ~ 2 L ff . 

The above calculations do not include pressure gradient effects. The pressure gradient 
terms has an all but negligible influence on the thin oil film solutions presented. According 
to equation 3.18, the pressure gradient associated with a saddle characterized by r x = ax 
and t z = bz , will have a magnitude of P x = — (3a+fc)/ tan#, where 6 ~ 30 deg is the angle of 
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departure from the surface of the streamline originating from the saddle into the flowfield. 
An order-of-magnitude analysis, shows that the oil film peak will be displaced upstream 
of the saddle by L p ~ ( 3 4- b/z)h$. For figure 13, the oil film peak associated with the 
saddle of separation would be displaced about L v % 4ho « 20 microns. This displacement 
caused by the pressure gradient would be about 1 grid point upstream. For the saddle of 
attachment shown in figure 13, the oil film peak would be displaced upstream even less, by 
L p % 10 microns, or 1/2 grid point. In the absence of pressure gradient terms, the saddle 
problem is symmetric in x as well as in 2 about the saddle. The surface tension terms are 
easily incorporated for the symmetric problem. However, the pressure gradient terms lead 
to a nonsymmetric (in x ) saddle problem. Because the effect of pressure gradient on the 
oil film is so small for these thin oil films, the pressure gradient terms were not included 
in the calculations for the oil film at this time. 

In the remaining portion of this saddle flow subsection, we use the above direct solu- 
tions in lieu of experimental data to explore several techniques for analysis of experimental 
oil film data. The context of the present study suggests three approaches. The first and 
simplest experimental thin oil film analysis in the vicinity of a saddle is to examine the 
surface streamlines for a rough estimate of b/a , thereby determining whether the saddle 
is a saddle of separation or a saddle of attachment. A second analysis is to examine the 
oil film thickness, h vs x/xo, along the centerline, z = 0, for estimates of the ratio of the 
shear gradients, b/a, and, again, whether the saddle is a saddle of separation or a saddle 
of attachment. However, the third and most thorough experimental thin oil film analysis 
of a saddle is to measure the thin oil film distribution h(x,z ), at one or more times, along 
with a surface streamline measurement, 7 (x, z), and then use an inverse solver, such as 
described earlier in this paper, to deduce the complete wall shear field, (r x (x, z), r z (x, z)), 
in the vicinity of the saddle. 

Placement of a pattern of oil dots or other conventional surface oil flow visualization 
techniques in the vicinity of a saddle point is straightforward and will result in an image 
similar to figure 3a or 3b. From figure 3a or 3b, for each selected surface streamline, s. 
we may identify an arbitrary point (x s , z s ) through which that streamline passes. Other 
arbitrary points ( Xi,Zj ) which lay on the same streamline may then be tabulated. In the 
vicinity of the saddle, the streamline coordinates will behave according to equation 3.19, 
rewritten here in logarithmic form as: 


log(— ) = (b/a) log(— ) (5.16) 

Z S X s 

A log-log plot of ( zjz s ) vs (x/x s ) for the saddle cases solved above is shown in figure 
14. Note that all the streamlines for each saddle case fall on one line associated with 
the b/a ratio for that case. The several saddle cases are clearly distinguishable in figure 
14, providing estimates of b/a and identifying which case is a saddle of separation or 
attachment. On a plot such as figure 14, each quadrant is associated with either saddles of 
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separation or with saddles of attachment but not both. The simplicity of this streamline 
analysis approach is compelling. Anticipated experimental complications not evident in 
this brief example include the fact that in general the sepratrix emanating from the saddle 
will typically be curved rather than straight. 

Figure 15 explores the utility of a plot of h vs logio(x/xo), as mentioned in the section 
on analytical solutions (Section 3). In figure 15, x is the distance from the saddle along 
the plane of symmetry, xq is the distance to the leading edge of the applied oil, and h is 
the oil film thickness. 

The ability to determine experimentally whether a given saddle is a saddle of sepa- 
ration or a saddle of attachment may prove desirable. We also may wish to estimate the 
ratio of the shear stress gradients, b/a. Figure 15 demonstrates how this might be done. 
Consider, in particular, the line shown in figure 15, which represents the special case where 
a -(- b = 0. This is a saddle which is intermediate between a saddle of separation and a 
saddle of attachment. For the case where a 4- b = 0, a plot of h vs logio(x/xo) will be a 
straight line, as indicated by the analytical solution given by equation 3.23 repeated here: 

ht = H = (— n 0 /a) ln(xo/x), for x < xo, a < 0, and a + b = 0 (5-17) 

For a saddle of separation, a + b < 0, a plot of h vs logio(x/x 0 ) will curve upward from 
the straight line as shown in figure 15, even with the effect of surface tension included. 

For a saddle of attachment, a + b > 0, a plot of h vs logio(x/xo) will curve downward 
from the straight line as shown in figure 15, again even with the effect of surface tension 
included. Clearly from figure 15, to use a log plot to unambiguously distinguish between a 
saddle of separation and a saddle of attachment, oil film thickness data must be obtained 
quite close to the saddle (x < 0.03xo). 

Inverse Solver, With and Without Surface Tension Terms 

An important aim of the present work was the development and validation of a 2D thin 
oil solver suitable for determining the 2D wall shear stress field on the surface bounding a 
complex 3D flow. To test, in this subsection, the inverse solver methods described earlier 
in this paper, we use the oil film thickness distributions obtained above by our direct 
solver for the 3D saddle of separation as the input to the inverse solver. A first oil film 
height distribution input case is considered where surface tension effects are not included 
in the direct solver, and then a second oil film height distribution input case where surface 
tension effects are included in the direct solver. Both the Two-Time-Level (eq. 4.26) and 
the One-Time-Level (eq. 4.30) versions of inverse solver algorithm are considered. Surface 
tension terms may also be explicitly added to the inverse Two-Time-Level solver algorithm 
and are also considered. 

Figure 16a and 16b show the T x (N/m 2 ) = — 2000x(ra) and the r z (N/m 2 ) = 1500z(ra) 
distributions, respectively, for the saddle of separation problem being considered. The 
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oil film distribution for this case was obtained using the direct solver without the surface 
tension terms (also see fig. 12) and was then used as the input to both the One-Time-Level 
and the Two-Time-Level (without surface tension) inverse solvers. For the Two-Time-Level 
solver, oil film distributions at t = 75 and t = 100 seconds were used, while for the One- 
Time-Level inverse solver, only the oil film distribution at t = 100 seconds was used. The 
wall shear stress distributions obtained by the inverse solvers are also shown in figures 16a 
and 16b. Note that uneven spacing in 2 was used to demonstrate the ability of both the 
direct and inverse solvers to deal with uneven grid spacing and also so as to better resolve 
the immediate vicinity of the saddle. Figure 16c shows the r x vs x results in the vicinity of 
the saddle to a finer scale than that of figure 16a. From the earlier discussion of the direct 
solver, we know the oil film distribution in the absence of surface tension agrees quite well 
with the analytical solution for this case. Figures 16a-16c clearly indicate that the inverse 
solver, both the One-Time-Level and the Two-Time-Level algorithms, also performs nearly 
ideally, in the absence of surface tension, in calculation of the 2D wall shear stress field 
about a saddle of separation. 

Surface tension effects do occur about a saddle of separation in an actual experimental 
setting. Thus, to examine the ability of the present inverse solver to deal with surface 
tension effects, we also consider the case where the input oil film thickness distribution 
was produced using the direct solver which included surface tension effects (see also the 
discussion of fig. 13). The ability of the direct solver to correctly incorporate surface 
tension effects was previously validated in the discussion of the surface tension bubble 
problem earlier in this section. 

Figures 17a and 17b are wall shear stress results, r x vs x and r 2 vs x. respectively, 
for the inverse solvers without surface tension terms in the inverse algorithm. Unlike 
the results in figures 16a-16c, however, the input oil film distribution for figures 17a and 
17b does include surface tension effects. Some error in the wall shear stress t x appears 
in the immediate vicinity of the saddle of separation as a result of introducing surface 
tension effects into the oil film thickness distribution; however, errors in t z close to the 
saddle become quite large. The surface-tension-induced errors are limited to a region of 
Ax ~ 2 L a about the saddle. In spite of the errors we observe, both the One-Time-Level 
and the Two-Time-Level inverse algorithms provide nearly identical results. 

Next, we add surface tension terms to the inverse Two-Time- Level solver to determine 
if these surface-tension-induced errors can be removed. Since the One-Time-Level inverse 
solver algorithm relies exclusively on the time self-similar relation h H/t, which is not 
valid for the surface tension term, the One-Time-Level inverse algorithm cannot be easily 
corrected for surface tension. 

Figure 18a and 18b are the wall shear stress results, t x vs x and r z vs x , respectively, 
for the Two-Time-Level inverse solver with surface tension terms added. Note that the 
surface-tension-induced errors close to the saddle of separation are still present but are 
smaller than those of figures 17a and 17b. To improve the accuracy for large time steps, an 
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assumption that h ~ \/t was incorporated into the derivation of both the direct solver and 
the inverse solver algorithms. Where the wall shear stress terms dominate, this assumption 
aids accuracy considerably. However, where surface tension terms become dominant, large 
time steps are not as accurately accomplished with this assumption. The direct solver 
results are obtained with a time step of At = 0.1 second. The inverse solver algorithm 
is designed for use with experimental data having large time steps between images and is 
tested here with a time step of At = 25 seconds. 

We have in this subsection demonstrated that the direct solver provides an excellent 
representation of the oil film thickness distribution about a saddle of separation or a saddle 
of attachment even to the extent of including surface tension effects in the vicinity of the 
saddle. Surface tension is found to influence only a quite limited region in the immediate 
vicinity of the saddle (Ax = 2 L a ~ 200 microns). The One-Time-Level and Two-Time- 
Level inverse solver algorithms have also been demonstrated to provide accurate wall shear 
stress distributions for the saddle of separation case, except for this limited surface-tension- 
dominated region. 
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6. CONCLUDING REMARKS 


The flow of a thin film of oil placed on a surface and then subjected to the 2D wall 
shear stress generated by a 3D aerodynamic flow is extensively considered. The governing 
partial differential equation is derived to include the effect of a wall shear stress and wall 
pressure gradient on the thin film of oil. Surface tension and gravity terms are included. 
We refer to the resultant partial differential equation as the “Thin Oil Film Equation.” 
The “Squire form” and the “Tanner form” of the thin oil film equation are shown, in 
Appendix A, to be equivalent through a metric transformation. Numerous analytical time 
self-similar solutions for the thin oil film equation are described. Of particular interest is 
the flow of a thin oil film in the vicinity of 3D aerodynamic surface streamline singularities 
such as the saddle of separation and the saddle of attachment. 

Another result of this study is that both direct and inverse numerical solution tech- 
niques are developed. In addition to the wall shear stress terms, these direct and inverse 
solution techniques may additionally include the effects of the wall pressure gradient, grav- 
ity and surface tension terms. A Two- Time-Level inverse solver is described which includes 
these effects. A One-Time-Level version of the inverse solver is for use in the absence of 
the wall pressure gradient and surface tension terms. When the oil film thickness varia- 
tions are provided from experimental images of the oil film, the inverse methods provide 
a rigorous mathematical basis for an improved form of the oil film based wall-shear-stress 
instrument. 

The numerical solvers are applied to several model problems having known analytical 
solutions so as to evaluate the fundamental accuracy of these numerical methods. Both the 
direct and inverse numerical solvers are shown to be stable, accurate and computationally 
efficient. One alternate tv sx solution method based on the local slope of the film thickness 
sometimes suggested elsewhere for oil film based wall-shear-stress instrumentation is shown 
to have fundamental accuracy problems leading to, in one example considered, errors of 
50% or more. A model surface tension problem with an analytical solution is also derived 
to demonstrate the ability of the present solvers to correctly implement the surface tension 
terms, which can become important near separation. 

An extensive application of the direct and inverse numerical solvers to the case where 
oil film flows on the surface in the vicinity of a saddle of separation provides a clear demon- 
stration of the success and utility of these numerical methods. The saddle of attachment 
case is also successfully treated. 

Techniques for the rigorous analysis of the behavior of thin oil films on test surfaces 
bounding complex 3D aerodynamic flows have been demonstrated. Instrumentation for 
the measurement of the 2D wall shear stress which use these numerical techniques should 
be accurate and suitable for use in complex 3D aerodynamic flows. 
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APPENDIX A 

The Equivalence of the Tanner Form and the Squire Form 
of the Thin Oil Film Equation 


In this appendix, we show that the 2D thin oil film equation of the form as given 
by Squire (see eq. 2.2) and of the form as given by Tanner (eq. A.6 derived below) 
are equivalent and may be derived from each other by means of a metric transformation. 
Understanding of the details of this transformation is useful in the construction of boundary 
conditions. The Tanner form of the thin oil film equation can be useful in formulation of 
certain types of boundary conditions, whereas the Squire form of the thin oil film equation 
is more compatible with current computational fluid dynamics numerical procedures. 

Consider that the test surface covered by an oil film may be described by a cartesian 
x-z coordinate system or by an s-n coordinate system attached to the surface streamlines 
as depicted in figure Al. Thus, n is constant for a given surface streamline, and s is the 
distance along the surface streamline. The local transformation between these two coordi- 
nate systems is accomplished by first a rotation ((x, z ) => (x, z)) by the local shear stress 
angle, 7, followed by a stretching dn = adz to account for the divergence (or convergence) 
of the surface streamlines. Note the stretching function, a, varies with position on the test 
surface. 

Clearly, 


ds ds . 

ds = —dx + — dz =cos7 dx + sin jdz 
ox dz 

, dn , dn . . , 

dn = —ax + — dz = — a sm'ydx + a cos 7 dz 
dx dz 


Inversion of the transformation matrix above gives: 


dx dx 

dx = —ds -I- —dn =cos7 ds — (sin 7/0 t)dn 
ds dn 

dz — —ds + —dn =sin7<2s + (cos'y/a)dn 
ds dn 


Further, transformation of partial derivatives between these two coordinate systems 
may be locally accomplished by: 


d ds d dn d 

dx dx ds dx dn 

d ds d dn d 

dz dz ds dz dn 


and 
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d _dx d dz d 

ds ds dx + ds dz 
d _dx d dz d 

dn dn dx + dn dz 

In the x-z coordinate system, consider the thin oil film equation of the form given by 
Squire (restricted to wall shear stress effects): 


dh d ,r x h 2 d ,r z h 2 

dt dx v 2 fj, a ; dz v 2 Ho 


)=0 


(A.1) 


To transform to the s-n coordinate system for the form given by Tanner, we first 
expand: 


d f r x h 2 d ,r z h 2 dh 2 /2/j, 0 . dh? j2yi 0 

— (it— ) + -St(t — ) =t cos 7 ^ — + rsin 7 - 


dx 2n 0 dz 2 Ho 


dx 


dz 


2 s/dr cos7 9rsin7 s 

+ {h 2 /2yi 0 ){ — — T - + — ^) 


dx 


dz 


Next, apply the coordinate transformation: 


d r x h 2 d r z h 2 dh 2 /2/x 0 . dh 2 /2fi 0 

— (77—) + — (tt— ) =t cos 7 (cos 7 — asin 7 — — — ) 


dx v 2 fi c 


dz 2ji 0 


ds ' dn 

, - , . dh 2 /2fi 0 dh 2 /2fj, 0 

+ r sm 7 (sin 7 h a cos 7 


+ (/i 2 /2^ g )(cos7 


9s 
dr cos 7 

ds 


) 


asm 7 


dn 
dr cos 7 


dn 

dr sin 7 9 rsin 7 , 

+ sm 7 — — (- a cos 7 — 77 — ) 


ds 


dn 


Now, contract the left hand side (mostly, cos 2 -I- sin 2 = 1 ), and rearrange giving: 


d r x h 2 d r z h 2 d ,rh 2 . l2/ 

= r(T-) + (rh /2 iXo) 


dx ' 2 ^ dl ^’ = a? 2 ^ ^ K ™ '^ 0)U dL 

which then leads to the s-n form for the 2 D thin oil film equation: 


dh drh 2 j2\i 0 
dt ds 


4 - (rh 2 /2fj. 0 )a^- = 0 
on 


(A.2) 


Note that in analyzing surface streamline images, a would seldom actually be estab- 
lished since the stretching, a, distorts angles, and the rotated coordinate, z, would be used 
locally rather than n. Thus, the substitution into the above equation of d'y/dz = ad'y/dn 
would be appropriate, particularly in the formulation of numerical boundary conditions: 
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dh drh 2 / 2/x 0 
dt + ds 


+ (rh 2 / 2^o) §J = 0 


(A. 3) 


To proceed toward Tanner’s form of the 2D thin oil film equation, consider figure A2 
which seeks to relate the streamline divergence term used by Tanner (note Tanner uses 
the symbol n rather than 77) to the shear stress angle 7 used in this paper. In figure A2, 
two surface streamlines are initially separated by a small distance, r) n = a(s,n)r]. Note, 
T] n is given in term of the stretched coordinate, n, whereas 77 is given in terms of the local 
unstretched coordinate, z. Thus, the angle formed between the two streamlines is given by 
A7 = (d'y/dz)rj. A small distance, As, along the streamlines the separation will increase 
by A77 = (d'y/dz)t]As. In the limit of small 77 and As, we obtain the relationship: 


1 <977 d'y _ c?7 

77 ds dz a dn 


(AA) 


Substituting this relation into equation A. 2 we obtain: 



drh 2 y ' 2 fx 0 

ds 


+ (t/i 2 /2// 0 )(— — 0 
77 ds 


Rearranging, we obtain Tanner’s form for the 2D thin oil film equation: 


(A. 5) 


dh 1 dryrh 2 / 2 fi 0 

dt 77 ds 


(A.6) 


Note, 77 can be obtained by integrating equation A.4 along a surface streamline (the 
choice for 770 at sq can be arbitrarily small): 


r, = Vo e { fso d ^ dids) (A. 7) 

Thus, the two streamlines initially some small distance apart, 770, at sq will be separated 
by 77 at s. But, in terms of the n coordinate, the separation, rj ni is fixed: 


which leads to: 


Tj n = q(s, n)rj(s,n) — constant 


a = a 0 e- { fso d i /dids) (A. 8 ) 

Equation A. 8 does give us the means to analyze a surface streamline image, with 
known 7 (x, z), for the stretching function, a(x, z), and known rotation function, 7(2:, z), 
required to establish the s-n coordinate system if so desired. 

The demonstration presented here that Tanner’s form and Squire’s form of the 2D thin 
oil film differ only in the coordinate system reveals the details of the required coordinate 
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system transformation. These details should prove useful in correctly formulating general 
boundary conditions. 


53 



APPENDIX B 
Numerical Forms for d 3 h/dx 3 

A centered numerical approximation to the partial derivative, d 3 h/dx 3 , may be found 
on an unevenly spaced grid ( Xi- 2 iXi-i,Xi,Xi + \,Xi + 2 ) by: 

( d 3 h/dx 3 )i — a.ihi +2 4- A^i+i 4- tihi 4- 4- 2 

4~ Xj_ 1 4- Xj_2 3Xj) 

Q/ • --I-. ' ' ' — — - ■ - ■ n -.-i - -- 

{x t + 2 ~ Xi){Xi+2 ~ Xi+i)(x i+2 ~ Xi-l)(Xi+ 2 “ 2) 

^ 6(Xj_j_2 4" Xi — 1 4” — 2 

(^t + l 2-i-bl ) (^i+l Xi— 1) Xi— 2) 

6(^+2 ^ x i-2 ~ 3Xi) 

(*^t — 1 Xi)(%i+2 l)(-^t— 1 “^*+1 )(^z— 1 — 2) 

^ 6(^+2 + x%+\ 4- Xj-\ — 3xj) 

(Xi-2 ~ X{)(Xi+2 ~ Xi- 2 )(Xi- 2 ~ Xi+ i)(Xi- 2 ~ ^i-l) 

ti = — (ai 4- Pi 4- 7i + ^i) 

A biased numerical approximation to the partial derivative d^h/dx 3 , may be found 
on an unevenly spaced grid (xi-\,Xi*Xi+\,Xi+ 2 ,Xi+z) by: 

(d 3 h/ dx 3 )i = 02^2+3 4- Pihi +2 4" 7t^i+i 4“ 4" 

— 6(Xi+2 4" #i+l 4" #1— 1 — 3x*) 

Q — — — 

(24+3 - Xi)(xi + 3 - Xi +2 )(a:i + 3 - 24+1X24+3 - Xj_i) 

q 6(24+3 4- Jj+1 4- Xj_i — 3Xj) 

(24+2 - X t )(x l+3 - X i+2 )(Xj + 2 - X i+1 )(x i+ 2 ~ Xi- 1) 

6(24+3 4- 24+2 4- Xj_ i — 3 Xi) 

= ~ 

3 *Ei+l ) (*Ei+l l) 

^ ___ 6(X2-t-3 4- X2+2 ^i+l ~ 3^2) 

( Xi — 1 Xi^){^Xi-\s Xi — X ) (X{ — 1 X{^2^)(Xi — 1 *£i+l) 

e 2 == — (ai 4- pi + 7i + 5») 
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Figure 2a. Saddle of attachment surface and plane of symmetry flowfield streamlines. 




Figure 2b. Saddle of separation surface and plane of symmetry flowfield streamlines. 
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Figure 4. Box-Implicit solution molecule. 




Figure 5. Thin oil film leading edge. 




Figure 6. Control volume for Finite- Volume algorithm. 
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Figure 7a. Oil film thickness from Box-Implicit direct solver 
(symbols) compared to analytical solution (lines) for 
linear wall shear stress (= 20 + 100 x) case. 
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Figure 7b. Oil film thickness from Box-Implicit Flux-Limit direct 
solver (symbols) compared to analytical solution (lines) 
for linear wall shear stress (= 20 + 100 x) case. 
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Figure 7c. Oil film thickness from Finite-Volume Upwind-Implicit 
direct solver (symbols) compared to analytical solution 
(lines) for linear wall shear stress (= 20 + 100 x) case. 
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Figure 7d. Oil film thickness from Finite-Volume Upwind-Implicit 
Flux-Limit direct solver (symbols) compared to 
analytical solution (lines) for linear wall shear stress 
(= 20 + 100 x) case. 
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Figure 7e. Oil film thickness from McCormack direct solver 
(symbols) compared to analytical solution (lines) for 
linear wall shear stress (= 20 + 100 x) case. 
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Figure 8a. Wall shear stress from Box-Implicit inverse solver 
Analytical oil film thickness at t=80 and 100 
seconds used as input to inverse solver. 
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Figure 8b. Wall shear stress from Box-Implicit inverse solver. 

Oil film thickness at t=80 and 100 seconds from 
Box-Implicit direct solver used as input to inverse solver. 
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Figure 8c. Wall shear stress from One-Time-Level Box-Implicit 
inverse solver. Box-Implicit direct solver oil film 
thickness at t=40 and 100 seconds as inverse solver input. 
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Figure 8d. Wall shear stress from local-slope method. 

Analytical oil film thickness at t=100 seconds used 
as input to local-slope method. 
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ickness at discrete times 
0 seconds) for constant wall 
onic initial thickness problem. 
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Figure 9d. Finite- Volume Upwind-Implicit direct solver oil film 
thickness at discrete times for constant wall shear 
stress, non-monotonic initial thickness problem. 
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Figure 9e. Finite- Volume Upwind-Implicit, Flux-Limit direct solver 
oil film thickness at discrete times for constant wall 
shear stress, non-monotonic initial thickness problem. 
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Figure 10. Initial and final oil film thickness, surface tension problem. 
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- - Analytical Solution, no surface tension 
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Direct Solver, no surface tension 
Direct Solver with surface tension 
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Figure 12. Oil film thickness for saddle problem, Box-Implicit 
method, t= 100 seconds. 
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Figure 13. Oil film thickness for detailed region near saddle 
point of saddle problem, Box-Implicit method, 
t=100 seconds. 
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Figure 14. Log quadrant analysis of saddle streamlines. 

Square symbols form figure 3a, a saddle of attachment. 
Circle symbols form figure 3b, a saddle of separation. 
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Figure 15. Centerline log analysis for saddle oil film thickness. 
□ Saddle of separation without surface tension terms. 
O Saddle of separation with surface tension terms. 

■ Saddle of attachment without surface tension terms. 
♦ Saddle of attachment with surface tension terms. 



Figure 16a. Oil film wall shear stress for saddle problem using 
Two- and One-Time-Level inverse solvers, 
no surface tension. 





Figure 16 b. Oil film wall shear stress for saddle problem using 
Two-Time-Level inverse solver, no surface tension. 
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Figure 16c. Oil film wall shear stress for detailed region near 
saddle point using Two- and One-Time-Level 
inverse solvers, no surface tension. 
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Figure 17a. Wall shear stress for saddle problem with surface 
tension effects using Two- and One-Time-Level 
inverse solver without surface tension terms. 
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Figure 17b. Wall shear stress for saddle problem with surface 
tension effects using Two- and One-Time-Level 
inverse solver without surface tension terms. 

Open symbols are Two-Time-Level inverse solver results. 
Filled symbols are One-Time-Level inverse solver results. 
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Figure 18a. Wall shear stress for saddle problem with 
surface tension effects using Two-Time-Level 
inverse solver with surface tension terms. 
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Figure 18b. Wall shear stress for saddle problem with 
surface tension effects using Two-Time-Level 
inverse solver with surface tension terms. 
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Figure A 1 . Streamline coordinate system. 
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Figure A2. Streamline divergence. 
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